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1.  INTRODUCTION 


Understanding  the  atomistic  mechanisms,  energetics  and  dynamics 
underlying  the  interactions  and  physical  processes  which  occur  when  two 
materials  are  brought  together  (or  separated)  is  fundamentally  important  to 
basic  and  applied  problems  such  as  adhesion,  contact  formation,  surface 
deformations,  materials  elastic  and  plastic  response  characteristics, 
materials  hardness,  microindentation,  friction,  wear  and  fracture.  These 
considerations  have  motivated  for  over  a  century  extensive  theoretical  and 
experimental  research  endeavors  of  the  above  phenomena  and  their  technological 
consequences.  Most  theoretical  approaches  to  these  problems,  with  a  few 
exceptions,  have  been  anchored  in  continuum  elasticity  and  contact  mechanics. 
Similarly,  until  quite  recently  experimental  observations  and  measurements  of 
surface  forces  and  the  consequent  materials  response  to  such  interactions  have 
been  macroscopic  in  nature. 

The  recent  emergence  and  proliferation  of  proximal  probes,  in  particular 
tip-based  microscopies,  sensitive  to  nano-  and  subnano-  meter  scale  structures 
provides  compelling  opportunities  for  studies  of  these  structures  which  are 
key  to  the  science  base  of  many  venerable  technological  problems.  In  addition 
these  probes,  coupled  with  advances  in  theoretical  understanding  of  the 
energetics  and  interaction  mechanisms  in  materials  and  the  development  of 
computer-based  materials  modeling  and  simulation  techniques,  open  new  avenues 
for  exploration  of  new  scientific  concepts  and  novel  materials  properties  and 
processes  on  the  subnanometer  scale. 

Our  research  program  is  aimed  at  developing  the  methodology  and 
theoretical  models  and  implementing  computational  techniques  to  allow 
investigations  of  the  aforementioned  materials  phenomena  on  the  atomistic 
level.  To  obtain  insights  into  the  energetics  and  dynamics  of  such  complex 
systems  and  phenomena  we  developed  molecular  dynamics  programs,  specifically 
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designed  for  studies  of  the  properties  of  materials  under  stress  and  for 
investigations  of  the  physical  consequences  of  tip-to-substrate  interactions. 

The  work  reported  herein  represents  new  avenues  for  studies  in  the  area 
of  nanomechanics  and  nanodynamics  as  well  as  for  investigations  using  the  new 
tip-based  spectroscopies.  As  we  describe  below  (under  "Summary  of  Research 
Achievements"),  in  addition  to  the  theoretical  studies,  we  have  initiated 
collaboration  with  experimental  groups  and  thus  have  developed  a  direct 
confrontation  between  theory  and  experiments. 
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2.  SUMMARY  OF  RESEARCH  ACHIEVEMENTS 


The  research  program  supported  under  the  above  AFOSR  grant  is  a 
continuation  of  a  research  effort  which  was  initially  supported  by  DARPA,  in 
the  context  of  an  initiative  directed  toward  investigations  of  tribological 
phenomena.  Our  studies  focus  on  basic  understanding  on  the  atomic  level  of 
tribological  phenomena  and  Of  the  properties  and  response  of  materials  under 
stressed  conditions.  In  order  to  gain  such  knowledge,  on  the  microscopic 
scale,  we  have  developed,  implemented  and  analyzed  computer-based  large-scale 
molecular  dynamics  simulation  methods,  which  allow  probing  of  the 
aforementioned  phenomena  with  refined  spatial  and  temporal  resolution. 

In  addition  to  the  theoretical  studies,  we  have  initiated  and  formed 
active  collaborations  with  experimental  research  groups  (in  particular  at  NRL 
and  Hughes)  which  result  in  a  comprehensive  understanding  of  fundamental 
processes  which  are  key  to  further  studies  and  to  future  technological 
applications. 

The  research  acheivements  of  our  investigations  may  be  summarized  as 
follows: 

(i)  Development  and  implementation  of  the  methodology  and  numerical 
simulation  techniques  for  studies  of  stress,  strain  and  finite  deformations  in 
materials.  See:  M.  W.  Ribarsky  and  Uzi  Landman,  Phys.  Rev.  B^8,  9522  (1988). 

(ii)  Development  of  the  methodology  and  simulation  methods  for  studies  of  the 
structural  and  dynamical  consequences  of  the  interaction  between  material  tips 
and  solid  surfaces.  These  studies  are  of  fundamental  importance  for 
understanding  the  atomic-scale  mechanisms  underlying  tribological,  adhesion, 
wear  and  other  interfacial  phenomena,  as  well  as  for  the  development,  design 
and  interpretation  of  experiments  using  the  emerging  tip-based  microscopies 
(STM,  and  in  particular,  Atomic-Force  Microscopy  ( AFM) ) . 

In  our  initial  studies  (see  Uzi  Landman,  W.  D.  Luedtke  and  M.  W. 
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Ribarsky,  J.  Vac.  Sci.  Technol.  £7,  2829  (1989);  reprint  attached)  we  have 
investigated  tip-substrate  interactions  for  covalently  bonded  materials 
(silicon),  demonstrating  the  local  nature  of  the  deformation  mechanisms,  the 
atomic  origins  of  the  stick-slip  phenomenon,  and  the  formation  of  connective 
necks  upon  separation  of  the  tip  from  the  substrate.  In  these  studies, 
simulations  of  atomically  sharp  tips,  blunted  ordered  and  disordered  tips,  in 
both  the  constant-height  and  constant-force  scanning  modes,  have  been 
performed. 

Our  recent  studies  (see  Uzi  Landman,  W.  D.  Luedtke,  N.  A.  Burnham  and 
R.  J.  Colton,  Science,  April  27  (1990);  preprint  attached)  focus  on  the 
atomistic  mechanisms  and  dynamics  of  adhesion,  nanoindentation  and  fracture  in 
metallic  systems.  To  simulate  faithfully  the  properties  and  response  in 
metals  one  must  employ  interaction  potentials  which  are  many-body  in  nature, 
including  the  electronic  contribution  to  the  cohesiveness  of  these  materials. 
In  our  studies  we  have  used  the  embedded  atom  theory  (EAM)  potentials  where 
the  cohesive  interactions  are  represented  by  embedding  functions,  supplemented 
by  parametrized  pair-interactions  which  describe  intercore  repulsion.  The 
particular  system  which  we  have  investigated  is  a  large  nickel  tip  (Ni, 
exposing  a  (001)  facet)  interacting  with  a  gold  (Au  (001))  substrate,  chosen 
because  of  the  difference  between  the  elastic  and  mechanical  properties  (such 
as  elastic  moduli  and  hardness)  of  these  materials. 

In  order  to  avoid  finite  size  effects,  we  have  used  in  our  molecular 
dynamics  studies  a  very  large  number  of  atoms  which  resulted  in  one  of  the 
most  challenging  simulations  ever  performed  (see  editorial  in  Science,  27 
October  (1989),  p.  4M5). 

The  results  of  our  studies  (see  full  details  in  the  attached  Science 
preprint)  reveal  the  atomistic  mechanisms  which  underiy  the  formation,  upon 
tip-to-substrate  approach,  of  an  intermetallic  adhesive  Junction,  monolayer 
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fracture  upon  tip-to-substrate  separation,  plastic  deformation  of  the  metal 
surface  region  upon  slight  indentation,  and  the  formation  and  elongation 
processes  of  atomically  thin  connective  necks  upon  separation  following 
contact. 

The  theoretical  predictions  from  our  studies  have  been  verified  by 
atomic-force  microscopy  experiments  where  the  predicted  hysteresis  in  the 
force  versus  tlp-to-sample  distance  relationship  was  measured.  In  addition  to 
guiding  the  interpretation  of  experimental  data  the  theoretical  results 
provide  significant  insights  into  the  microscopic  processes  of  deformation, 
adhesive  friction,  fracture  and  wear  phenomena,  as  well  as  opening  new  avenues 
for  critical  assessment  and  reformulation  of  traditional  continuum-based 
theories  of  materials  deformation  and  response  in  the  elastic,  plastic  and 
inelastic  regimes. 

(iii)  Another  facet  of  our  research  program  focuses  on  the  properties  of 
ionic  solids  under  stress  and  on  interfacial  phenomena  in  ionic  materials. 

In  particular,  the  system  which  we  have  investigated  is  crystalline  CaF2, 
which  for  the  bulk  material  exhibits  the  onset  of  superionic  conductivity  at 
~1420K  (while  the  melting  point  is  at  ~1700K).  Due  to  this  property  we  expect 
that  the  mechanical  properties  of  this  material  will  show  pronounced  variation 
with  temperature  (as  is  indeed  observed  experimentally).  Our  simulations  of 
the  response  of  bulk  crystalline  CaF2  to  applied  strain,  for  a  wide  range  of 
temperatures,  show  a  strong  dependence  of  the  resolved  shear  yield  stress  on 
temperature. 

Since  our  main  interest  is  in  interfacial  systems  we  have  conducted 
extensive  investigations  of  the  surfaces  of  this  material,  exposing  the  (110) 
or  the  (111)  crystalline  faces.  Our  investigations  show,  for  the  first  time, 
the  onset  of  surface  superionic  conductivity  at  temperatures  much  below  that 
of  the  bulk.  Furthermore,  the  onset  temperature  for  surface  superionic 


5 


conductivity  (Ts.  )  depends  strongly  on  the  exposed  crystalline  face.  Thus 
while  for  the  (111)  surface  of  CaF2  only  a  slight  enhancement,  over  bulk 
values,  of  mobility  of  the  F"  anions  is  observed  at  temperatures  below  the 
onset  of  bulk  superionic  conductivity  (t|\  ~  1400K) ,  the  (110)  surface 

becomes  a  superionic  conductor  at  a  temperature  as  low  as  1000K  (i.e.,  400K 
below  the  bulk);  see  Figure  at  the  end  of  this  section. 

In  addition  to  the  above  findings  we  observed  no  premelting  of  the 
surfaces  ((111)  and  (110))  at  temperatures  below  the  bulk  melting  temperature. 
Thus  the  Ca+2  sublattice  maintains  it's  crystalline  nature  up  to  the  melting 
transition  of  the  crystal. 

It  is  of  interest  to  note  in  this  context  that  the  region  of  space  above 
the  topmost  layer  of  the  surface  does  not  serve  as  a  diffusion  channel. 
Instead,  the  presence  of  the  surface  serves  to  lower  the  activation  energy 
barrier  for  generation  of  anion  vacancies  which  in  turn  promotes  jump- 
diffusion  processes  in  the  top  layer  of  the  substrate  and  in  deeper  layers 
at  higher  temperatures.  Finally  the  magnitude  of  superionic  conductivity  (as 
measured  for  example  by  the  F~  jump  rate)  at  the  surface  of  the  material 
(below  and  above  the  bulk  onset,  Tgic  ~  1400K)  is  noticably  larger  than  in 
the  bulk. 

These  results  indicate- that  this  system  offers  a  potentially  most 
interesting  material,  in  which  tribological  characteristics  are  coupled  to  a 
unique  physical  process  (i.e.,  super  ionic  conductivity).  Indeed  our  current 
studies  of  the  yielding  processes  of  the  material  at  elevated  temperatures, 
and  simulations  of  calcium  difluoride  tips  interacting  and  scanning  over  a 
CaF2  crystalline  surface  show  a  very  pronounced  dependence  on  the  substrate 
temperature  as  well  as  on  the  rate  of  the  scans.  In  a  way  one  may  postulate 
at  this  point,  that  the  interface  is  capable  of  "self- lubricating”  due  to  the 
enhanced  mobility  of  F  at  the  interface.  Further  investigations  on  this 
system  are  in  progress  and  will  be  reported  when  completed. 


6 


F  JUMP  RATE  (psec -1) 


BULK  e*4 


+  '  V 


1000.0 1100.0 1200.01300.01400.01500.01600.01700.' 

TEMPERATURE  Ck) 


p 

51  ^  ^  OV\sd~  hulk  5u^l|evii'c 

^A+O.  Co«duc^i+fc  1  J,  .1,  ,  /, 


S  I  C 


■jov'  ]~\SL  OV\sdr  hulk  5uf£Jt  I'evii'c 


—  Fi<p.  2-Cl  — 


F  JUMP  RATE  (psec -1) 


5uv-fac>£  /aj&'t  (111) 


« 

O 


T 

O 


CN 

o 


o 

o 

1000.0 1100.0 1200.01300.0M00.01500.01600.01700.0 

TEMPERATURE  (k) 


Fitf,  zb  — 


F  JUMP  RATE  (psec.-D 


i  i  \  i  i  i  i  i 

1000.0 1100.0 1200.01300.01400.01500.01600.01700.C 


TEMPERATURE  (k) 


SU'T'fdi 

Note  iU  of  c owduc'frvrf^  Jxz/fr«j~  Tstc 

(C&MpcflKG  ts-  Ftg>2-*) 

-  F/pc  - 


F  JUMP  RATE  (psec-1) 


Surface  jVjfirfllO) 


TEMPERATURE  (k) 


F  JUMP  RATE  (psec -1) 


SfryfflCl.  yH  (110) 


o  o 


o  o  ° 


1000.0 1100.0 1200.01300.01400.01500.01600.01700.C 

TEMPERATURE  (k) 


-  fty  ie  - 


3.  CONFERENCES  AND  SYMPOSIA  PRESENTATIONS 

1.  Invited  talk  (U.L.)  at  the  ONR  Workshop  on  "New  Theoretical  and 
Computational  Approaches  to  the  Structure  and  Properties  of  Cracks  in 
Solids",  NIST,  Washington,  D.C.,  September  1989. 

2.  Invited  talk  (U.L.)  at  International  Workshop  on  "Many-Atom  Interactions 
in  Solids",  Pajulahti,  Finland,  June  1989. 

3.  AFOSR  Molecular  Dynamics  Contractors  Meeting,  Captiva  Island,  Florida, 
November  1989. 

4.  Invited  talk  (U.L.)  at  the  Gordon  Conference  on  "Fundamentals  of 
Tribology",  New  Hampshire,  June  1990. 

5.  Invited  talk  (U.L.)  at  the  "International  Conference  of  Scanning  Micro¬ 
scopy  and  Nanoscience",  Baltimore,  Maryland,  July  1990. 
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pieces  that  will  leave  a  residue. 

Having  dismissed  poivols,  Jung  turned  to 
a  new  base  stock  whose  composition  is 
pioprietary  information.  Although  this  ma¬ 
terial  does  decompose  at  700°F,  it  leaves 
little  residue.  Further,  although  the  oil 
bums,  it  bums  so  slowly  that  oil  loss  is  not  a 
problem,  says  John  Fairbanks  of  the  Depart¬ 
ment  of  Energy,  who  oversees  the  joint 
Akzo/Cummins  program.  Only  the  thin  film 
of  oil  along  the  top  ring  of  the  piston  at  the 
top  of  its  stroke  bums,  so  the  oil  losses  are 
not  much  more  than  the  losses  of  a  conven¬ 
tional  oil  in  a  conventional  engine. 

Along  with  the'  base  stock,  Jung  also 
developed  a  package  of  additives  based  on 
theoretical  considerations.  For  instance,  he 
says,  inorganic  chemistry  can  be  used  to  pre¬ 
dict  which  products  arc  going  to  attach  them¬ 
selves  most  firmly  to  the  metal  engine  parts 
and  thus  serve  as  useful  antiwear  additives. 

Cummins  has  already  tested  the  synthetic 
lubricant  that  Jung  developed,  and  so  far  it 
looks  like  it  will  not  only  be  a  good  high- 
temperature  oil  but  that  at  lower  tempera¬ 
tures  it  will  be  far  superior  to  conventional 
oils,  Fairbanks  says.  In  one  test,  Cummins 
ran  the  oil  for  300  hours  at  a  top  tempera¬ 
ture  of  850°F,  and  it  had  virtually  no  ash 
deposit  and  relatively  little  oil  loss,  Fair¬ 
banks  says.  Jung  adds  that  it  could  be  avail¬ 
able  commercially  within  a  few  years  and 
that  it  will  be  priced  about  the  same  as 
current  synthetic  automobile  oils,  such  as 
Mobil- 1. 

As  a  bonus,  the  new  lubricant  should  help 
reduce  automotive  pollution,  Fairbanks 
says.  “In  conventional  engines,  30%  to  60% 
of  the  particulate  emissions  come  from  lu¬ 
bricants,”  he  says.  But  the  Akzo  lubricant 
“should  volatize  direedy  into  a  gas  and  not 
form  particulates,  which  will  play  a  major 
role  in  reducing  emissions.” 

At  the  Midwest  Research  Institute,  Sutor 
also  decided  that  hotter  engines  would  de¬ 
mand  a  new  approach  to  designing  lubri¬ 
cants.  “We  got  tired  of  trying  what  you  can 
get  off  the  shelf,”  Sutor  says.  So  2  years  ago 
he  went  “back  to  first  principles”  and  de¬ 
signed  a  lubricant  “from  a  dean  sheet  of 
paper.” 

Like  Jung,  Sutor  decided  additives  were 
to  blame  for  much  of  the  deposit  formation 
and  designed  a  new  set  of  additives.  Unlike 
Jung,  he  decided  to  try  to  make  his  base 
stock  more  resistant  to  oxidation.  For  in¬ 
stance,  since  longer  chain  lengths  make  the 
molecule  more  prone  to  oxidation,  he  chose 
a  base  stock  with  shorter  chain  lengths. 

Then  he  went  to  chemistry  to  see  if  he 
could  protea  the  sites  on  the  base  stock 
molecules  that  were  particularly  prone  to 
oxidation.  The  oxidative  attack  comes  at  a 
particular  spot  on  the  chain,  he  discovered. 
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For  most  of  its  history,  the  “science  of  rubbing  bodies”  has  been  an  empirical 
profession — more  engineering  than  predictive  science.  “There  is  no  unified  theory  of 
tribology,”  savs  Irwin  Singer,  a  tribologist  at  the  Naval  Research  Laboratory.  And 
although  tribologists  believe  they  understand  friction,  wear,  and  lubrication  qualita¬ 
tively,  they  cannot  usually  make  numerical  predictions  except  by  relying  on  models 
derived  from  experimental  data. 

In  the  past  few  years,  several  theorists  have  set  out  to  remedy  this  situation.  Their 
goal  is  to  understand  enough  about  what  happens  at  a  molecular  level  when  two 
surfaces  move  against  each  other  to  be  able  to  suggest  new  directions  and  new 
approaches  on  the  experimental  side  of  tribology.  And,  although  the  models  they  are 
developing  cannot  yet  handle  the  complex  interactions  between  surfaces  and  lubri¬ 
cants,  the  researchers  believe  they  are  laying  the  groundwork  for  true  science. 

For  example,  Uzi  Landman  and  David  Luedtke.  two  physicists  at  Georgia  Tech, 
have  used  molecular  modeling  on  a  Cray  supercomputer  to  predict  what  will  happen 
when  the  tip  of  a  thin  nickel  needle  is  pushed  down  onto  a  flat  gold  surface  and  then 
lifted  up  again.  In  their  model,  the  needle  tip  contained  2,000  nickel  atoms,  while  the 
surface  consisted  of  8,000  gold  atoms.  TTie  goal  was  to  calculate  the  individual 
trajeaories  of  all  10,000  atoms,  which  Landman  and  Luedtke  did  in  two  steps. 

They  first  used  basic  quantum  mechanics  to  determine  the  electrostatic  potential  in 
which  each  atom  moves.  “All  of  the  questions  of  adhesion,  cohesion,  making  of 
bonds,  breaking  of  bonds — these  are  really  quantum  mechanics  questions,”  Landman 
says.  Once  the  researchers  knew  the  potential,  they  could  calculate  the  motion  of  each 
atom  by  a  straightforward  application  of  Newtonian  physics.  Straightforward,  but 
definitely  not  simple.  With  10,000  atoms  involved,  the  supercomputer  takes  up  to 
150  hours  to  calculate  a  motion  that  in  real  time  takes  only  a  billionth  of  a  second. 

The  computer  simulation  shows  a  quite  peculiar  action.  As  the  nickel  tip 
approaches  the  gold  surface,  the  attractive  force  on  the  tip  increases  slowly.  Then 
suddenly,  when  the  tip  is  about  4  angstroms  from  the  surface,  the  top  layer  of  gold 
atoms  on  the  surface  breaks  loose  and  bulges  out  to  touch  the  tip.  The  gold  atoms  in 
the  top  layer  of  the  surface  are  more  strongly  attracted  to  the  nickel  atoms  in  the  tip 
than  to  their  fellow  gold  atoms.  Landman  plains.  If  the  tip  is  now  retraacd,  it  pulls 
a  thin  neck  of  gold  atoms  along  with  it.  mter  this  neck  breaks,  there  is  a  one- atom- 

thick  coating  of  gold 
atoms  on  the  nickel 
tip  and  the  gold  sur¬ 
face  is  left  with  one 
very  thin  column  of 
atoms  pointing  up. 

Nancy  Burnham 
and  Richard  Colton 
at  the  Naval  Research 
Laboratory  subse¬ 
quently  used  an  atom¬ 
ic  force  microscope  to 
check  these  calcula¬ 
tions  by  observing  the 
interaction  between  a 
nickle  needle  and  a 
gold  surface  directly. 
And,  Landman  says, 
the  behavior  they  measured  was  exactly  as  predicted  by  the  supercomputer. 

Such  simple  theoretical  models  may  be  useful  in  understanding  such  chings  as  how 
fractures  form,  Landman  says,  but  the  greater  value  will  appear  once  researchers 
model  much  more  complicated  systems.  “The  stage  is  set  for  a  foil  interaction  between 
theory  and  experiment,”  Landman  says.  Eventually,  he  hopes  that  a  theoretical 
understanding  of  the  molecular  forces  between  two  surfaces  or  between  a  surface  and 
a  lubricant  will  give  scientists  a  better  understanding  of  what  is  going  on  at  the 
macroscopic  level.  "To  make  progress  in  the  molecular  engineering  of  lubricants,”  he 
says,  “you  need  to  know  the  molecular  details  of  the  process.”  ■  R.P. 


Atomic  attractions.  As  the  nickel  tip  pulls  away  jrom  the  gold 
surface,  it  takes  a  layer  of  gold  atoms  (lighter  spheres)  with  it. 
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Formulations  and  methodologies  of  molecular-dynamics  simulations  of  material  systems  evolving 
under  applied  finite  external  perturbations  are  developed  and  discussed.  Focusing  on  interfacial  sys¬ 
tems,  composed  of  interfacing  crystalline  solids  characterized  by  differing  interatomic  interactions 
and  atomic  sizes,  the  mechanisms  and  dynamics  of  response  and  stress  relief  in  the  elastic,  plastic, 
and  inelastic  regimes  are  investigated.  Critical  values  of  the  external  perturbations  (stress  and 
strain)  are  determined,  showing  dependence  on  the  nature  of  the  interface  and  the  ambient  condi¬ 
tions  (thermally  adiabatic  versus  isothermal). 


L  INTRODUCTION 

Basic  understanding  of  the  structure  and  dynamics  of 
materials  and  their  properties  often  requires  knowledge 
on  a  microscopic  level  of  the  underlying  energetics  and 
interaction  mechanisms  whose  consequences  we  observe 
and  measure.  The  structural,  mechanical,  and  dynamical 
response  of  material  systems  to  externa]  stresses,  on  the 
atomic  level,  are  key  issues  in  developing  a  fundamental 
understanding  of  a  number  of  systems  and  phenomena  of 
coupled  basic  and  technological  interest  such  as  tribolo¬ 
gy,  lubrication,  wear,  material  fatigue  and  yield,  crack 
propagation,  stress-induced  phase  and  structural  trans¬ 
formations,  and  hydrodynamical  phenomena,  to  name 
just  a  few.  Although  the  above-listed  phenomena 
represent  everyday  experiences  and  have  been  observed 
and  studied  for  a  long  time,  detailed  microscopic  theories 
of  them  (with  few  exceptions)  are  lacking.  Nevertheless 
large  bodies  of  empirical  data  and  in  some  cases  phenom¬ 
enological  model  descriptions  have  been  developed. 

Common  observations  related  to  the  above  phenome¬ 
na,  and  of  thermomechanical  properties  and  response  of 
materials  in  general,  are  usually  made  at  the  continuum 
level.1  Consequently  (and  naturally)  the  development  of 
theoretical  understanding  of  these  phenomena  followed 
the  "continuum  modeling”  approach.1  The  methodology 
of  the  development  of  these  models  is  based  on  the  princi¬ 
ples  of  mass,  momentum,  and  energy  balance,  and  the 
formulation  of  constitutive  equations.  While  the 
mathematical  formulation  of  classical  models  of  the 
mechanical  response  of  matter  (such  as  the  classical 
theories  of  elasticity  and  hydrodynamics)  achieved  a  high 
degree  of  sophistication,  current  focus  is  on  the  incor¬ 
poration  of  the  knowledge  about  the  microscopic  behav¬ 
ior  of  materials  in  continuum  models  which  attempt  to 
describe  macroscopic  observations.  This  is  done  via  the 
introduction,  into  the  continuum  models,  of  a  set  of  state 
variables  which  provide  averaged  (coarsened)  representa¬ 
tions  of  the  relevant  microscopic  quantities.  In  addition 
we  should  note  that  most  applied  models  of  mechanical 
response  are  limited  to  the  elastic  (small  spatial  deforma¬ 
tion)  and  linear  (small  rates  of  application  of  the  external 


perturbation)  regimes.  Attempts  to  incorporate  inelastic 
response  and  nonlinear  effects  result  in  a  great  (often 
prohibitive)  complexity. 

Computer  molecular-dynamics  (MD)  simulations2-9 
where  the  evolution  of  a  physical  system  is  simulated, 
with  refined  temporal  and  spatial  resolution,  via  a  direct 
numerical  solution  to  the  model  equations  of  motion  are 
in  a  sense  computer  experiments  which  open  new  avenues 
in  investigations  of  the  microscopic  origins  of  material 
phenomena.  These  methods  alleviate  certain  of  the  major 
difficulties  which  hamper  other  theoretical  approaches, 
particularly  for  c  molex  systems  such  as  those  character¬ 
ized  by  a  large  number  of  degrees  of  freedom,  lack  of 
symmetry,  nonlinearities,  and  complicated  interactions. 
In  addition  to  comparisons  with  experimental  data,  com¬ 
puter  simulations  can  be  used  as  a  source  of  physical  in¬ 
formation  which  is  not  accessible  to  laboratory  experi¬ 
ments,  and  in  some  instances  the  computer  experiment  it¬ 
self  serves  as  a  testing  ground  for  theoretical  concepts. 

At  this  stage,  the  development  and  application  of  com¬ 
puter  simulations  to  studies  of  materials  phenomena,  and 
in  particular  the  mechanical,  structural,  and  dynamical 
response  of  materials  to  external  perturbations,  can  be 
used  in  order  to  investigate  the  microscopic  mechanisms 
of  material  response  and  to  determine  the  critical  materi¬ 
al  parameters  (cohesive  energies,  interaction  potentials, 
atomic  sizes,  crystallographic  structure)  and  ambient 
conditions  (for  example,  thermally  adiabatic  versus  iso¬ 
thermal  conditions  and  dependencies  on  the  strength  and 
rate  of  applied  external  perturbations)  which  govern  the 
observed  trends.  The  information  which  can  be  obtained 
via  microscopic  simulations  could  then  guide  the  formu¬ 
lation  of  continuum  models,  provide  numerical  values  for 
certain  physical  material  characteristic  parameters  which 
could  be  incorporated  in  such  models,  and  provide  gui¬ 
dance  for  the  design  and  analysis  of  well-controlled  ex¬ 
perimental  investigations. 

Traditionally,  MD  simulations  have  been  employed  in 
studies  of  systems  of  fixed  shape  and  size  of  the  periodi¬ 
cally  replicated  calculational  cell  (i.e.,  constant  volume 
simulations).  More  recently  methods  for  simulating  sys¬ 
tems  in  which  the  volume  and  shape  of  the  calculational 
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cell  may  vary  dynamically  have  been  developed,10-15 
opening  the  way  to  investigations  of  a  large  number  of 
materials  phenomena  in  which  the  dynamical  freedom  of 
the  system  to  change  volume  and/or  structure  (or  phase) 
is  essential.15  In  addition,  a  number  of  methods  have 
been  developed  for  simulations  of  flow  and  hydrodynami- 
cal  systems,16-19  which  allow  detailed  investigations  of 
these  nonequilibrium  phenomena.  Using  molecular- 
dynamics  techniques  various  studies  of  the  mechanical 
properties  of  solids  and  fluids  have  been  reported. 
Among  the  investigations  of  solid  systems  we  note  studies 
of  stressed  solids  and  crack  propagation,20-25  dislocation 
energetics  and  dynamics,23  structural  transformations  in 
crystal  lattices  under  uniaxial  tension  or  compres¬ 
sion,1 1,15,26,27  sliding  and  migration  of  grain  boundaries,15 
simulations  of  plastic  deformations  and  shock  wave  dy¬ 
namics,26,29  studies  of  stressed  interfaces,17  and  calcula¬ 
tions  of  elastic  coefficients  using  MD  simulations.30 

In  this  paper  we  investigate,  using  MD  simulations  the 
microscopic  dynamical  response,  deformation,  and 
stress-relief  mechanisms  at  crystalline  solid  interfaces 
subject  to  externally  applied  perturbations.  Studies  of 
solid  interfaces  are  of  inherent  as  well  as  applied  interest 
in  the  development  of  atomistic  models  of  friction,  solid 
lubrication,  and  wear  phenomena,  since  plastic  deforma¬ 
tion  which  occurs  in  the  vicinity  of  the  interface  between 
sliding  materials  is  the  principal  mechanism  for  the  dissi¬ 
pation  of  frictional  work.1  The  objectives  of  these  studies 
are  (i)  identification  of  the  mechanisms  for  solid  interfa¬ 
cial  systems  of  deformation,  stress  accumulation  and  re¬ 
lief,  and  the  dynamical  response  to  external  perturba¬ 
tions,  (ii)  identification  of  the  dependence  of  the  above 
phenomena  on  material  characteristics,  such  as  bonding 
strength,  atomic  sizes,  and  interface  crystallography,  and 
on  ambient  conditions  (thermally  adiabatic  versus  iso¬ 
thermal),  and  (iii)  the  development  and  critical  assess¬ 
ment  of  MD  simulation  methods  for  investigations  of  the 
above  phenomena.  In  Sec.  II  we  discuss  molecular- 
dynamics  methodologies  and  techniques  for  simulations 
of  finite  deformations  of  material  systems  under  stress. 
The  setup  of  the  systems  and  results  of  simulations  using 
MD  simulation  methods,  under  several  thermally  adia¬ 
batic  and  isothermal  conditions  and  for  systems  of 
different  sizes  and  material  characteristics  are  described 
and  compared  in  Sec.  III.  A  summary  of  our  findings  is 
given  in  Sec.  IV. 

II.  MOLECULAR-DYNAMICS  FORMULATION 

In  the  molecular-dynamics  (MU)  method  the  equations 
of  motion  for  a  set  of  interacting  particles  are  integrated 
numerically  and  properties  of  the  system  are  obtained 
from  the  generated  phase-space  trajectories  of  the  sys¬ 
tem.2-9  In  the  case  of  an  extended  system  periodic 
boundary  conditions  (PBC’s)  are  imposed  while  in  simu¬ 
lations  of  finite  aggregates6  no  such  device  is  used.  The 
starting  point  of  a  MD  simulation  is  a  well-defined  micro¬ 
scopic  description  of  the  physical  system,  in  terms  of  a 
Hamiltonian  or  a  Lagrangian  from  which  the  equations 
of  motion  are  derived.  In  the  early  applications  of  the 
MD  method  to  extended  systems  the  simulations  were 


performed  in  the  microcanonical  ( E,V,N )  ensemble,  i.e., 
the  energy  £,  volume  V,  and  number  of  particles  A  in  the 
calculational  cell  (which  is  periodically  repeated)  are  con¬ 
stants.  The  desire  to  investigate  physical  situations  in 
which  the  volume  and  shape  (geometry)  of  the  material 
system  may  vary  (in  response  to  an  external  pressure  or 
applied  stress)  as  well  as  circumstances  where  isothermal 
(i.e.,  constant  temperature  T)  pertain,  motivated  the  de¬ 
velopment  of  new  MD  formulations,  such  as  the  ansatz 
Lagrangian,10-13  constrained  dynamics,4,9  and  the  Nose- 
dynamics  methods.9  In  the  following  we  review  first  the 
ansatz-Lagrangian  method,  as  extended  by  Parrinello  and 
Rahman11  (PR)  for  calculations  of  systems  under  exter¬ 
nally  applied  stresses,  and  its  extensions  to  treat  finite  de¬ 
formations.12  Subsequently  we  discuss  an  alternative 
method  for  dynamical  studies  of  deformations  in  which 
the  system  evolves  under  an  imposed  constant-strain  rate. 
Finally,  we  describe  a  variant  of  the  MD  method  for 
finite  deformation  studies  and  a  general  method  which  is 
particularly  suited  for  investigations2,18,19  of  systems  in 
which  flow  phenomena  may  occur. 

A  The  ansatz-Lagrangian  method:  The  ensemble 

In  the  PR  formulation11  the  particles’  coordinate  vec¬ 
tors  whose  components  are  ria  (where  /  =  1,2 . A  is  a 

particle  index  and  a  =  1,2,3  is  a  coordinate  index)  are  ex¬ 
pressed  in  terms  of  scaled  coordinate  vectors  sia  (where 
Ogsia<  1,  /  =  1,2, . . .  ,N  and  a= 1,2,3)  and  a  matrix  U 
whose  columns  H,,  H2,  and  H3  are  the  vectors  A,B,C 
which  span  the  edges  of  the  MD  calculational  cell, 

r,«=  2  Bafit  <*'  . ^  a,/?=  1,2,3) .  (1) 

0-i 

The  volume  of  the  cell.  Cl,  is  given  by  det(#  )=  |  |  H  |  | . 

Regarding  the  components  of  the  H  matrix  as  dynami¬ 
cal  variables,  Parrinello  and  Rahman11  have  originally 
proposed  a  Lagrangian  (a  Hamiltonian  formulation  can 
also  be  developed12,27)  from  which  the  equations  govern¬ 
ing  the  time  evolution  of  the  3 A  +9  degrees  of  freedom 
(3A  particle  and  nine  calculational  cell  degrees  of  free¬ 
dom)  are  obtained: 

£pR=i  2  K* ft*, ))]+£«„- u*,, , 

«-i 

(2) 

where  dotted  variables  indicate  differentiation  with 
respect  to  time,  the  superscript  T  denotes  a  transpose,  the 
metric  tensor  Q  is  defined  as  Q=H  rfl,  U  is  the  potential 
energy  of  the  particles  (which  depends  on  the  specific 
form  of  the  interparticle  interaction),  and  the  cell  kinetic 
energy  was  chosen  by  PR  as 

=  (3) 

where  IP  is  a  mass  parameter.  For  the  case  of  a  system  to 
which  an  external  stress,  specified  by  the  stress  tensor  q_, 
is  applied  the  potential  energy,  t/Mn,  is  given  by 

l/cell  =  ft0Tr(2,£)  ,  (4) 
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where  the  strain  tensor  £  can  be  written  in  terms  of  the 
metric  tensor  Q  and  the  reference  value,  /£0,  of  /£,  in  the 
form 

£=|[(a0"1)rfiS-1-l]  •  (5) 

Using  Eq.  (5)  in  Eq.  (4),  and  denoting  by  Pe  =  j Tr<e^ )  the 
external  hydrostatic  pressure,  the  expression  for  in 

the  small-strain  limit"  [i.e.,  Tr(£)~(ft  —  n0)/ft0]  is 

C/^-P.lft-IV+ITriSfi),  (6) 

where 

I=n0£0“l(2*-^I><£“1)7'  (7) 

and  ft0=det(2i0).  Note  that  for  the  case  of  an  isotropic 

external  stress  £=0,  the  second  term  in  Eq.  (6)  vanishes 
and  only  the  contribution  from  the  hydrostatic  pressure, 
Pt,  remains. 

The  equations  of  motion  which  are  derived  in  the  usual 
manner  from  the  above  Lagrangian  [Eq.  (2)]  for  the  parti¬ 
cles  are 

N 

mfii  =  —  2  (i  =1,2, . . .  ,N)  , 

j-i 

(8a) 

where  Xij=ri~HdU/driJ),  and  those  describing  the 
dynamical  motion  of  the  calculational  cell  are  given  by 

W&M{L-P,\.)A-UZ  ,  (8b) 

where  A  is  the  “area  tensor,”  J4=n(fl-1)r  [i.e., 
,^5=30/3 Hap)],  and  the  elements  of  the  internal  mi¬ 
croscopic  stress  tensor  are  given  by 


N 

N 

2 

PiaPi0^mi~  2  Xtjrtjeftjfi 

i-i 

J-‘  +  • 

(a,0=l,2,3)  (8c) 

where  =  |  r;  —  ry  [  and  the  momentum  p( =mlflS(  . 

The  above  formulation  generates  an  isoenthalpic- 
isostress  (9i,ot,N)  ensemble  is  the  enthalpy)  with  the 
elastic  energy  given  by  [see  Eq.  (6)] 

Etl  =Pt(Cl  -n0)+ ft0Tr[(£,  -P,l)£] 

=P<(n-n0)  +  |Tr(i£E)  .  (9) 

The  sum  of  Etl  and  the  particle  kinetic  and  potential  en¬ 
ergies  is  a  constant  of  the  motion  in  this  formulation. 11,12 

Two  issues  raised  at  this  point  are  (i)  the  appropriate 
choice  of  the  reference  state  #011,12  [see  Eq.  (5)j,  and  (ii) 
the  nonuniqueness12  of  the  form  for  £cen  [Eq.  (3)].  Con¬ 
sidering  issue  (i)  first  it  is  noted  that  as  far  as  the 
definition  of  the  strain  tensor  £  [Eq.  (5)]  is  considered  the 
choice  of  Ho  is  arbitrary."  12  However,  the  choice  of  the 
reference  system  is  of  importance  when  considering  the 
elastic  energy,  particularly  when  finite  deformations  are 
involved.12  An  expression  for  the  virtual  work  6 E  per¬ 
formed  in  a  virtual  deformation  (i.e.,  the  elastic  energy 
stored  in  the  medium)  of  a  deformable  medium,  whose 
volume  in  the  deformed  state  is  ft,  was  derived  by  Mur- 


naghan,31 

5 E=Eel  =  JTtUKoU ~ 's^H -XTHl )6f]dft  ,  (10) 

with  the  strain  tensor  £  defined  by  Eq.  (5)  and  the  un¬ 
stressed  system,  H0,  as  the  reference  state.  The  integral  is 
to  be  performed  over  the  final  state,  i.e.,  stressed  body. 
The  “calculational  cell  tensors”  Ho  and  H  for  the  un¬ 
stressed  and  stressed  system,  respectively,  occur  in  Eq. 
(10)  since  a  position  r  in  the  deformed  system  relates  to 
the  corresponding  position  r0  in  the  undeformed  body  via 
the  Jacobian  H Hq\  i.e.,  t=HHqX*o-  In  the  limit  of 
the  infinitesimal  theory  of  elasticity  Eq.  (10)  becomes 

Ee/  =  f[Tr(ze8£)dn  .  (11) 

For  a  uniform  system  undergoing  a  finite  deformation  Eq. 
(10)  takes  the  form 

£r/  =  ftoTr(i£)  ,  (12) 

where  the  “thermodynamic  tension”  x»  is  defined31,32  by 

I=-^flo£r,H«(fl",)7flor-  (13) 

For  the  infinitesimal  case  we  obtain 

Js^ftoTrfu.f)  .  (14) 

We  emphasize  that  in  Eq.  (12)  the  reference  state,  fl0,  is 
the  unstressed  system.  In  Eq.  (14)  however,  the  reference 
state  is  a  system  under  stress  and  as  proposed  by  PR,11  in 
this  case  Ho  could  be  chosen  as  the  average  of  the 
stressed  system,  <fl>.  Since  the  thermodynamic  tension 
X  is  the  quantity  appearing  in  the  thermodynamic  expres¬ 
sions,31’"  e.g.,  the  enthalpy  +ft<>Trti£),  where  E 
is  the  particle  energy,  Ray  and  Rahman12  proposed  the 
isoenthalpic-isotension,  ensemble  where  x  is 

constant  [but  not  see  Eq.  (13)].  For  this  ensemble  the 
equations  of  motion  for  the  particles  are  as  given  by  Eq. 
(8a)  and  the  equation  for  H  [Eq.  (8b)]  is  replaced  by 

W&=<lA-ILL,  (15a) 

where 

•  E=ft0ffo',I<flo'',>7'-  (15W 

In  some  of  the  simulations  discussed  in  the  next  section 
we  have  used  the  ( )  ensemble.  In  the  following 
section  we  formulate  an  alternative  dynamical  simulation 
method  in  which  the  system  evolves  under  a  constant- 
strain  rate. 

B.  Constant-strain-rate  simulations 

In  the  constant-strain-rate  (CSR)  simulations  the  calcu¬ 
lational  cell  deforms  in  a  prescribed  way  at  a  prescribed 
constant  rate.  In  this  method  only  the  particle  equations 
of  motion  need  to  be  solved  since  the  evolution  of  the  cal¬ 
culational  box  is  given,  i.e.,  the  components  of  H.  do  not 
evolve  freely.  We  remark  that  since  for  a  given  deforma¬ 
tion  (i.e.,  prescribed  strain)  the  internal  stresses,  SL,  see 
Eq.  (8c),  which  the  system  develops  (in  equilibrium,  or 
steady  state  in  case  of  flow)  in  response  to  the  deforma- 
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tion  can  be  calculated,  the  CSR  and  (7f,r,N)  ensemble 
(see  Sec.  II  A)  simulations  can  be  operated  “self- 
consistently”  by  using  these  calculated  stresses  as  applied 
external  stresses  in  simulations  where  the  dynamical  evo¬ 
lution  of  the  calculational  cell  is  included11,12  (i.e.,  using 
the  ensemble,  or  the  isoexternal-stress  method 

described  in  Sec.  C).  The  formulation  of  the  CSR  method 
is  facilitated  by  defining  the  streaming  velocity,  u,  as 
u=fis  .  (16) 

Using  this  definition  the  strain-rate  tensor,  £,  is  given  by 

£=Vu =(ft~l)TftT .  (17) 

The  equations  of  motion  for  the  particles  can  be  now 
written  [compare  to  Eq.  (8a)]  as 

£=Vu=(fl-')rflr.  (17) 

The  equations  of  motion  for  the  particles  can  be  now 
written  [compare  to  Eq.  (8a)]  as 

2  rP* •  (18) 

7- 1 

where  the  peculiar  momentum  p,  is  defined  by 

p  j—m&ki  (19a) 

and 

p,=m,(ffi,+fi)i,).  (19b) 

Rearranging  terms  and  using  Eqs.  (19),  we  write  Eq.  (18) 
as 

P/  =  -  2  V/;-£p<  (20) 

y-i 

and 

r,=p,/m,+£rr,  .  (21) 

We  note  that  Eq.  (20)  is  the  same  as  the  “Doll’s  tensor” 
equation  of  motion  introduced  by  Evans  and  Hoover.16 

As  discussed  by  Ray  and  Rahman,12  in  the  ansatz  La- 
grangian  [Eq.  (2)],  terms  involving  H a,  are  omitted.  Em¬ 
ploying  a  Hamiltonian  formulation  these  authors  de¬ 
rived12  an  alternative  equation  of  motion  which  includes 
the  ftt,  terms  [see  Eq.  (2.12)  in  Ref.  12], 

m/HVi  =  -  2  X<jTij~ — 2mifli,  .  (22) 

7-1 

For  a  process  in  which  the  system  is  deformed  at  a  con¬ 
stant  strain  rate,  £  is  constant  for  a  Couette  geometry 
(i.e.,  fluid  sheared  between  two  parallel  plates  held  at  a 
fixed  distance  from  one  another,  see  Sec.  Ill  C),  the  term 
involving  ft  in  Eq.  (22)  vanishes  and  using  the  definitions 
of  £  [Eq.  (17)]  and  the  peculiar  momentum,  p,  [Eq. 
(19a)],  the  above  equation  of  motion  becomes 

N 

Pi  =  -  2  xijrij  —£TPi  •  (23) 

7-1 

We  remark  that  this  equation  involves  the  transpose  of  £ 
[compare  to  Eq.  (20)],  and  is  the  same  as  the  local-rest- 
frame  dynamics  (“Sllod”)  equation  of  motion  (transpose 
of  the  non-Newtonian  dynamics  “Doll’s  tensor”)  suggest¬ 
ed  and  used  by  Evans  and  Morris.16 


Thus  we  observe  that  for  constant-strain  rate  the  equa¬ 
tions  of  motion  generated  via  the  Parrinello  and  Rah¬ 
man,"  and  Ray  and  Rahman12  formulations  reduce  to 
those  used  in  the  nonequilibrium-molecular-dynamics 
(NEMD)  method.16  In  the  constant-strain-rate  (CSR) 
simulations  described  below  we  have  used  Eqs.  (23),  fol¬ 
lowing  NEMD  studies.  However,  we  remark  that  simu¬ 
lations  employing  Eqs.  (20)  yield  qualitatively  similar  re¬ 
sults. 

C.  Isoexternal-stress  formulation 


In  the  previous  two  sections  we  discussed  methods  for 
molecular-dynamics  simulations  of  system  which  under¬ 
go  finite  deformations,  i.e.,  the  isoenthalpic-isotension, 
iJf,T,N)  ensemble,  and  constant-strain-rate  (CSR) 
methods.  In  order  to  allow  simulations  of  systems  under 
constant  external  stress,  g^,  we  have  developed  a  new 
method  which  in  addition  avoids  certain  pathological  cir¬ 
cumstances  which  may  be  encountered  when  using  the 
previous  methods,  particularly  when  flow  develops. 
While  we  do  not  use  this  method  of  simulations  in  this 
paper,  we  include  a  brief  discussion  of  it  for  complete¬ 
ness. 

In  order  to  develop  an  isoextemal-stress  dynamical 
method,  we  consider  the  Lagrange  equations  of  motion 
which  include  nonconservative  forces,”  Qi( 


d_ 

dt 


3  L 

_3 L_ 

a*/a 

3?.o 

—  Qia  * 


(24) 


(a=l,2,3)(/  —  1,2 . N+  9)  , 


where  the  set  of  q,’s  and  4,  ’s  include  the  dynamical  de¬ 
grees  of  freedom  of  the  particles  and  of  the  calculational 
cell.  The  virtual  work,  8 E,  done  by  the  nonconservative 
forces  is  given  by 

8£=2Q,-6q,  •  (25) 

i 

Recalling  the  definition  of  the  area  tensor,  A  [sec  Eq. 
(8b)]  we  note  that  the  components  (  A,,  A2,  A3)  of  A  are 
the  areas  of  the  calculational  cell  times  their  respective 
normals  (i.e.,  A1=H2XH3,  etc.).  Therefore  due  to  a  vir¬ 
tual  deformation  of  the  calculational  cell  (in  which  the 
vectors  H„  are  changed  by  8Ha)  each  of  the  surfaces  Aff, 
is  displaced  by  8Ha.  The  amount  of  virtual  work  done  in 
displacing  each  of  the  surfaces  due  to  the  force  a Aa  is 
8ft  Tg_t  Aa  (for  a  =  1 ,2, 3 ).  For  a  general  virtual  displace¬ 
ment  8ft  the  virtual  work  done  is 

8£  =Tr(8ftT a,, A)  »  (26) 


which  can  be  shown  to  be  the  same  as  the  result  given  in 
Eq.  (10)  for  a  uniform  system.  Comparing  Eqs.  (25)  and 
(26)  the  Lagrange  equations  of  motion  [Eq.  (24)]  for  the 
calculational  cell  take  the  form 


d_  3  L 

^  3 


-■=!£— 2 af.«r<*Mrp  (a,)3=  1,2,3), 

atla0  y 


(27) 


where  we  have  denoted  the  q' s  and  q's  corresponding  to 
the  calculational  box  by  H  and  H  and  the  possible  depen- 
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dence  of  the  external  stress  tensor,  q^,  on  time  is  indicat¬ 
ed.  Using  the  Lagrangian  given  in  Eq.  (2),  but  without 
the  term  which  is  explicitly  accounted  for  in  Eq. 
(27),  equations  of  motion  for  the  calculational  cell  tensor, 
ft  can  be  derived  yielding, 

Wft  =  (g_-(^)A  ,  (28) 

in  place  of  Eq.  (8b),  or  (15a)  in  the  previous  formulations. 
The  particle  equations  of  motion  remain  unchanged  [Eq. 
(8a)].  Note  that  here  £,  is  constant  unlike  the  case  of  the 
(9i,r,N)  ensemble  equations  of  motion  and  that  the  ques¬ 
tion  of  reference  state  does  not  arise.  The  system  howev¬ 
er  is  nonconservative. 

We  turn  next  to  the  issue  of  the  nonuniqueness  of  the 
form  of  the  calculational  cell  kinetic  energy,  K ^  [see 
Eqs.  (2)  and  (3)].  This  question  has  been  addressed  before 
and  it  has  been  shown  that  any  form  for  which  is  a 
function  of  H  and  H  generates  the  isoenthalpic-isobaric 
ensemble  for  sufficiently  large  N.n  Furthermore,  note 
that  for  such  choices  the  particle  equations  of  motion 
[Eq.  (8a)]  are  independent  of  as  are  also  equilibrium 
equal-time  ensemble  averages.16,12 

It  is  natural  to  demand  that  the  equations  of  motion 
governing  the  dynamical  evolution  of  the  system  will  be 
invariant  under  transformations  connecting  equivalent 
computational  cells19  (which  may  differ  in  shape  but  are 
of  the  same  volume,  and  contain  equivalent  sets  of  parti¬ 
cles).  Since  flslHj.Hj.Hj),  defining  the  computational 
cell,  fully  incorporates  the  periodicity  of  the  system,  i.e., 
any  two  equivalent  points  in  the  periodically  replicated 
system  may  be  connected  by  a  vector  n,H,-Hi2H2 
+n3H3  where  n,,  n2,  and  n3  are  integers,  we  define19  the 
class  of  transformations  under  which  the  system  is 
translationally  invariant  as  T  transformations.  These 
transformations  have  the  property  that  all  the  elements  of 
the  (3x3)  transformation  matrix  X are  integers  and  that 
det(  T)—  1  (to  assure  conservation  of  volume). 

It  has  been  recently  shown19  that  previous  choices  for 
K^n  [such  as  the  Parrinello-Rahman  form,11  Eq.  (3),  and 
variants  thereof12]  lead  to  equations  of  motion  for  the 
calculational  cell  which  do  not  obey  the  above  invariance 
requirement.  Moreover,  it  has  been  shown  that  using 
*cdi  by 

K^\W^t(ftATAftT)  ,  (29) 

generates  equations  of  motion  invariant  under  T  transfor¬ 
mations  and  in  addition  leads  to  satisfaction  of  the  virial 
theorem  at  equilibrium  (where  it  is  appropriate,  i.e.,  iso¬ 
tropic  pressure  a. 1*  =  \P,V- 

In  addition  to  providing  a  more  complete  description 
by  virtue  of  obeying  a  natural  invariance  requirement, 
this  formulation  presents  a  technical  advantage  which  is 
of  importance  particularly  in  simulations  of  flow  sys¬ 
tems.1®  Consider,  for  example,  a  fluid  flowing  under  the 
influence  of  an  external  stress  under  steady-state  condi¬ 
tions.  For  concreteness  take  the  external  stress  tensor 
to  be  a  symmetric  tensor  £rafl=C(8a|5pj-l-8a38p,),  where 
C  is  an  arbitrary  constant.  Under  the  influence  of  the 
external  stress  the  calculational  cell  will  eventually  be¬ 


come  extremely  nonorthogonal  (in  fact  due  to  the  absence 
of  effective  resistance  to  flow,  such  skewed  calculational 
cell  may  develop  in  fluid  simulations  even  without  exter¬ 
nal  stresses).  Suppose  that  the  interparticle  interaction 
potentials  in  the  system  extend  to  a  finite  range  smaller 
than  the  linear  dimension  of  the  calculational  cell.  Then 
in  calculating  the  force  acting  on  particle  i  we  need  to  in¬ 
clude  contributions  from  all  particles  located  inside  the 
range  of  interaction  in  the  calculational  cell  and  from 
their  images  in  the  first  (so-called  minimum  image) 
periodic  replications  of  the  cell.  As  long  as  the  calcula¬ 
tional  cell  does  not  deviate  much  from  being  cubical  (i.e., 
orthogonal)  the  search  for  periodic  images  can  be  imple¬ 
mented  most  efficiently  by  finding  the  nearest  integer  to 
each  component  of  the  scaled  coordinates  [see  Eq.  (1)], 
tj  —  a,- .  However,  for  very  skewed  cells  there  is  no  such 
efficient  procedure  and  since  this  segment  of  the  calcula¬ 
tion  is  one  of  the  most  time  consuming  parts  of  the  simu¬ 
lation  this  poses  a  serious  impediment.  A  solution  to  this 
problem  is  to  transform,  in  the  course  of  the  simulation, 
the  skewed  cell  to  an  equivalent  one  which  is  more  or¬ 
thogonal,  as  often  as  necessary.  This  transformation  in¬ 
volves  only  the  small  number  of  degrees  of  freedom  asso¬ 
ciated  with  the  calculational  cell  (ft).  Of  course,  in  order 
to  afford  such  transformations  one  must  require  that  the 
equations  of  motion  governing  the  system  evolution  be 
invariant  under  these  transformations,  which  is  achieved 
via  the  choice  of  given  in  Eq.  (29).  Recent  simula¬ 
tions  in  our  laboratory1®  of  stressed  systems  beyond  the 
yield  point  and  of  sheared  fluids  support  the  usefulness  of 
the  new  method. 

IIL  RESULTS 

To  investigate  the  properties  and  the  dynamics,  ener¬ 
getics  and  response  of  a  solid  to  externally  imposed  forces 
we  have  performed  molecular-dynamics  simulations  on  a 
model  system  at  different  ambient  conditions.  To  illus¬ 
trate  the  different  simulation  modes  discussed  in  the  pre¬ 
vious  section  we  compare  results  obtained  for  systems  in 
which  the  calculational  cell  was  allowed  to  respond 
dynamically  to  an  external  thermodynamic  tension  (see 
Sec.  II  A)  with  those  obtained  for  a  calculational  cell 
which  is  deformed  (keeping  the  volume  unchanged)  at  a 
constant-strain  rate  (the  CSR  method  described  in  II B). 
The  two  simulation  modes  illustrate  various  aspects  of 
the  deformation  process  and  the  approach  of  the  system 
to  yield.  The  dynamical  cell  mode  of  simulation  may  be 
likened  to  the  situation  where  a  load  and  shear  stress  are 
applied  to  materials  in  frictional  contact,  resulting  in  pro¬ 
cesses  which  are  accompanied  by  volume  changes,  while 
the  constant-strain-rate  mode  corresponds  to  experiments 
in  which  a  material  pressed  between  two  parallel  moving 
plates  (slabs,  which  are  kept  at  constant  distance)  is 
strained.  As  the  stress  on  the  solid  system  is  increased,  or 
equivalently  as  the  strain  is  increased,  the  materials  de¬ 
form,  first  elastically  and  then  inelastically  culminating  in 
yield  when  a  disruption  of  the  material  occurs.  In  this 
paper  we  investigate  the  response  of  the  solid  as  it  ap¬ 
proaches  the  yield  point.  Studies  beyond  yield  require  a 
modification  of  the  simulation  procedure  (see  Sec.  II C) 
and  will  be  discussed  elsewhere.1® 
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A.  System  setup 

In  this  study  we  focus  on  the  mechanisms  and  dynam¬ 
ics  of  the  response  of  a  system  containing  an  interface  be¬ 
tween  two  materials  to  external  stresses.  The  model  sys¬ 
tem  which  we  employ  in  our  simulations  consists  of  N 
particles  (NA  of  type  A  and  NB  of  type  B,  NA  +NB=N) 
interacting  via  pairwise  6-12  Lennard-Jones  (U)  poten¬ 
tials 

V(r)= 

(30) 


12 

r 

(a,0)=A,B  , 


where  A  and  B  represent  two  types  of  materials.  The 
solid  is  set  up  initially  in  a  face-centered-cubic  crystalline- 
structure,  with  (111)  layers  (N /NL  particles  per  layer) 
with  the  z  axis  along  the  [111]  direction,34  and  three- 
dimensional  periodic  boundary  conditions  are  used  (see 
Fig.  1).  The  well-depth  parameter  (e^)  of  the  interaction 
potential  between  particles  in  layers  1  -NLA  is  taken  to  be 
twice  that  for  particles  in  layers  NLA  +  1-N  (i.e., 
eAA  =2eBB )  corresponding  to  a  soft,  solid,  lubricating 
material  (the  B  system)  pressed  between  hard-material 
slabs  (the  A  system).  In  order  to  isolate  the  dependence 
of  the  system  response  on  the  interaction  strength  (poten¬ 
tial  depth)  parameter,  e,  from  that  due  to  differences  in 
the  atomic  sizes  (and  thus  interatomic  distances)  associat¬ 
ed  with  the  parameter  cr  in  the  LJ  potential,  we  have  per¬ 
formed  first  simulations  in  which  eAA—  2eBB  and 
cr  AA  =aBB,  and  then  simulations  in  which  <rBB 
=  1.5 a AA.  The  interspecies  LJ  potential  parameters 
were  chosen  as  oAB—(<rAA  +oBB)/2  and 
eAB— AAeBB^n-  In  the  following  we  use  reduced 
units35  where  energy  and  temperature  are  expressed  in 
units  of  eAA,  length  in  units  of  a AA,  stress  in  units  of 
(fAA  /°\a  )» and  the  time  unit  (t.u.)  is  (mA/eAA  ),/2<rAA. 
In  the  integration  of  the  equations  of  motion  we  used  a 
time  step  At  =0.0075  t.u.  where 

t.u .=(mA/€AA)in(TAA  , 

which  results  in  energy  conservation  (to  six  significant 
figures)  for  extended  runs. 

In  most  of  our  simulations  we  have  used  systems  con¬ 
sisting  of  N  =  1260  particles  with  18  (111)  layers  [NL) 
and  ArLA=^VZJ=9  (unless  specified  differently  results  are 
for  this  size  systems).  However,  in  order  to  investigate 
and  assess  system  size  dependencies  we  have  also  carried 
out  comparative  simulations  for  systems  containing 
IV  =  1890  particles  with  NLA  =9  (i.e.,  NA  =630)  and 
Nlb  =  18  (i.e.,  Nb  =  1260).  Some  results  from  these  simu¬ 
lations  will  be  exhibited  when  the  size  dependence  is  dis¬ 
cussed. 


B.  Applied  thermodynamic-tension 
simulations:  (9i,r, N)  ensemble 

1.  Adiabatic  conditions 

In  all  our  simulations  we  equilibrate  first  the  initial  sys¬ 
tem  at  a  reduced  temperature  of  7=0. 11  (note  that  a 


pure  bulk  U  crystal  melts  at  7=:  0.7,  and  thus,  since 
eBB=eAA  /2,  the  bulk  melting  point  of  the  soft  material, 
B,  is  half  that  of  the  A  material).  Subsequently  we  apply 
to  the  system  a  load  along  the  [111]  (z)  direction  (normal 
to  the  interface,  as  shown  in  Fig.  1)  at  a  rate  of 
0.0025 (e/a3)/At  until  a  load  value  of  0.5  is  reached. 
Following  equilibration  under  this  load  we  apply  to  the 
system  a  thermodynamic  tension  r  in  the  [  llO]  direction 
(see  Fig.  1)  at  a  rate  of  fxl=0.00\25(e/ai)/At  and  fol¬ 
low  the  evolution  of  the  system  until  it  fails  [to  keep  the 
system  against  a  rigid  rotation  a  symmetric  thermo¬ 
dynamic  tension  tensor  r  is  applied  (i.e.,  tu  =  riI )]. 

To  investigate  the  dependence  on  the  thermal  ambient 
conditions  we  distinguish  between  simulations  where  the 
system  is  thermally  isolated  during  the  application  of  the 
shear  [discussed  in  (1)]  and  simulations  where  isothermal 
conditions  are  maintained  [discussed  in  (2)].  The  tem¬ 
poral  history  of  the  applied  tension  versus  time  for  the 
adiabatic  and  isothermal  conditions  are  shown  in  Figs. 
2(a)  and  2(b),  respectively.  As  the  system  evolves  under 
the  applied  shear  it  develops  internal  stresses  which  are 
calculated  using  the  positions  and  forces  on  the  particles 
using  Eq.  8(c).  For  a  rate  of  increase  of  the  applied  exter¬ 
nal  shear  which  is  much  smaller  than  the  characteristic 
relaxation  rates  of  the  material  the  system  would  evolve 
on  a  phase-space  trajectory  which  corresponds  to  an 
equilibrium  (or  quasiequilibrium)  path  for  which  the 
values  of  the  external  and  internal  stresses  are  equal 
along  the  system  evolution.  For  externally  applied  ther¬ 
modynamic  tension  above  a  certain  critical  value  [rr  and 
the  corresponding  calculated  critical  external  stress  aK, 
see  Eq.  (13)]  failure  of  the  system  will  occur,  evidenced  by 
a  drop  of  the  internal  stresses  to  zero  (i.e.,  stress  relief) 
and  an  unbounded  variation  of  the  volume  of  the  system. 
Since  the  rate  of  increase  of  the  applied  thermodynamic 
tension,  r,  indicated  by  the  solid  line  (slanted)  in  Fig.  2,  is 
not  slow  enough  to  allow  the  system  to  relax  at  all  times, 
the  critical  values  thus  obtained  serve  only  as  rough  esti¬ 
mates  (upper  bounds)  to  the  true,  quasistatic,  critical 
values.  Therefore,  to  obtain  a  more  accurate  estimate  of 
the  critical  stress  we  have  performed  simulations  at  con¬ 
stant  values  of  the  thermodynamic  tension  indicated  by 


FIG.  1.  A  schematic  of  the  calculations!  cell.  NLA  is  the 
number  of  layers  in  the  A  (hard)  material  and  NLI  the  number 
of  layers  in  the  B  (soft)  material.  The  interface  is  between  layers 
Nla  and  Nla  + 1.  The  directions  of  the  applied  load  and  shear 
stresses  are  indicated.  Three-dimensional  periodic  boundary 
conditions  are  employed  in  the  simulations. 
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FIG.  2.  Summary  of  the  thermally  adiabatic  (a)  and  iso¬ 
thermal  (b)  simulations.  The  applied  thermodynamic 

tension  r„  vs  time  (in  reduced  t.u.)  is  shown.  The  slanted  solid 
line  represents  the  rate  of  application  of  ra  (following  equilibra¬ 
tion  under  load).  Horizontal  lines  correspond  to  simulations  at 
constant  r„,  at  the  values  given  by  the  intersection  of  the  lines 
with  the  r  axis.  The  termination  of  a  horizontal  line  by  a  sym¬ 
bol  ( X )  represents  that  at  the  corresponding  time  the  system 
yielded  in  response  to  the  applied  perturbation.  Lines  which  do 
not  terminate  by  an  x  correspond  to  r  values  for  which  the  sys¬ 
tem  did  not  yield.  From  these  simulations  the  critical  values  for 
structural  transformation  and  eventual  yield  (given  in  Table  I) 
were  determined.  Bold  solid  lines  correspond  to  the  system 
after  slip  and  stacking-fault  formation. 

the  horizontal  solid  lines  in  Fig.  2  (in  these  simulations 
the  system  is  brought  to  the  desired  value  of  r  and  then 
evolves  under  that  constant  value  of  the  thermodynamic 
tension).  In  this  figure  the  absence  of  a  symbol  (x)  at  the 
end  of  a  horizontal  line  indicates  that  the  system  did  not 
yield  for  that  value  of  r,  while  a  termination  of  the  line  by 
x  ii.dicates  yield  at  the  corresponding  time.  We  note  that 
as  the  value  of  r  under  which  the  system  evolves  is  in¬ 
creased  the  time  required  for  the  system  to  yield  shortens. 
The  bold  solid  lines  in  Fig.  2  correspond  to  the  system 
after  slip  and  stacking-fault  formation.  From  these  simu¬ 
lations  we  obtain  a  value  of  rc .  „  =0.86±0.02  for  the  crit¬ 
ical  thermodynamic  tension  (corresponding  to  crKXI 
=0.95±0.03),  see  Table  1.  Inspection  of  the  real-space 
trajectories  of  the  particles  reveals  that  yield  involves 


TABLE  I.  Critical-yield  values  of  the  thermodynamic  ten¬ 
sion  (rc ),  external  stress  (aK ),  calculated  via  Eq.  (13),  and  inter¬ 
nal  stresses  (ac ),  calculated  via  Eq.  (8c),  obtained  from  adiabat¬ 
ic  and  isothermal  simulations  using  the  (9i,r,N)  ensemble.  The 
r„  a„,  and  a,  values  correspond  to  the  critical  values  of  the 
above  quantities  for  structural  transformations  prior  to  yield 
(stacking  fault  or  slip,  as  discussed  in  the  text)  in  the  ( !H,r,N ) 
ensemble  simulations.  Under  the  constant-strain  heading,  re¬ 
sults  are  given  for  the  critical  strain  (y,)  to  bring  about  an  in¬ 
elastic  structural  transformation,  in  constrain-strain  simula¬ 
tions,  and  the  corresponding  values  of  the  internal  stresses,  (a, ) 
for  thermally  adiabatic  simulations  of  18-  and  27-layer  systems, 
as  well  as  results  for  isothermal  constant-strain  simulations  of 
the  18-layer  system.  All  quantities  are  in  Lennard-Jones  re¬ 
duced  units. 


Adiabatic 

Isothermal 

0.8  ±0.02 

0.84±0.03 

0.78±0.01 

0.83±0.03 

0.78±0.01 

0.84±0.04 

0.86±0.02 

0.95±0.01 

°tc.n 

0.95±0.03 

0.99±0.005 

0.94±0.02 

1.00±0.02 

Constant  strain 

Adiabatic  (18-layer) 

Isothermal  (18-layer) 

0.063 ±0.002 

0.065±0.003 

0.975 

1.02 

Y,/a, 

0.065 

0.064 

Adiabatic  (27-layer) 

r*» 

0.076±0.002 

1.14 

Y,/<r, 

0.067 

interplanar  motion  of  the  (111)  atomic  layers  (relative  to 
one  another)  in  accordance  with  observations  that  the 
main  operative  slip  system  in  fee  crystals  consists  of  (1 1 1) 
planes  in  the  [1 10]  direction. 

For  a  certain  range  of  values  of  r  the  system  undergoes 
structural  transformations  which  do  not  result  in  total 
yield.  This  is  indicated  by  the  bold  solid  horizontal  lines 
in  Fig.  2  and  the  corresponding  values  are  given  in  Table 
I  under  T^a  (0.8±0.2  for  the  adiabatic  system  and 
0. 84±0.03  for  the  isothermal  one).  To  investigate  the 
mechanisms  of  response  and  stress  relief  prior  to  yield  we 
have  performed  extensive  studies  of  the  behavior  of  the 
adiabatic  system  under  an  applied  ru=0.81.  It  should 
be  noted  that  removing  the  external  perturbation  prior  to 
the  onset  of  the  structural  transformation  results  in  the 
system  returning  to  the  average  unstressed  conditions, 
while  when  doing  so  past  the  structural  transformation 
the  system  remains  in  the  deformed  state.  Records  of  the 
average  temperature  (D,  potential  energy  (Ep ),  kinetic 
energy  ( Ek ),  and  elastic  energy  (Etl)  [the  work  done  on 
the  system,  see  Eq.  (12)]  versus  time  (starting  at  t  =  140 
t.u.,  the  intersection  of  the  horizontal  line  corresponding 
to  r=0.81  with  the  solid  line  in  Fig.  2)  up  to  the  yield 
point,  are  shown  in  Fig.  3.  [Note  that  the  sum  of  the  en¬ 
ergies  in  Figs.  3(b)— 3(d)  is  constant  ]  Inspection  of  the 
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FIG.  3.  Records  of  the  temperature  ( T ),  potential  energy 
(£, ),  kinetic  energy  (Ek ),  and  elastic  energy  (£w )  vs  time  (start¬ 
ing  at  t>140  t.u.  obtained  via  an  ( 3i,r,N )  simulation  at 
r„=0.81  under  thermally  adiabatic  conditions.  The  sharp 
changes  in  all  quantities  correspond  to  structural  transforma¬ 
tions.  All  quantities  are  in  Lennard-Jones  r<.ut— w  units. 


figures  reveals  that  the  average  temperature  of  the  system 
increases  from  the  initial  value  (in  reduced  LJ  units)  of 
0.11  achieving  a  new  value  of  —0.25  at  about  175  t.u. 
Correspondingly,  the  magnitude  of  the  system  average 
potential  energy  decreases  and  the  stored  elastic  energy 
increase  in  the  above  time  interval.  In  addition  we  find 
that  the  volume  of  the  system  increases  during  that  time 
period.  Inspection  of  the  structure  of  the  system  shows 
that  during  the  time  interval  —150-175  t.u.  the  system 
undergoes  structural  transformations,  occuring  in  the  re¬ 
gion  occupied  by  the  soft  (solid  lubricant)  material  (layers 
10-18).  The  positions  of  atoms  in  a  central  (112)  slice  of 
the  system  before  and  after  the  transformation  displayed 
in  Figs.  4(a)  and  4(b),  respectively,  show  clearly  the  for¬ 
mation  of  stacking-fault  region  [i.e.,  change  from  the  ab- 
cabc.  .  .  registry  sequence,  see  Fig.  4(a)],  in  the  soft  ma¬ 
terial  (open  circles).  In  addition  to  the  formation  of  the 
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FIG.  4.  Atomic  trajectories  for  a  (112)  central  slice  of  the 
system  simulated  with  r=0. 81:  (a)  at  (  =  139.1  t.u.  and  (b)  at 
(  =206.6  t.u.,  i.e.,  past  the  structural  transformation  (see  Fig.  3). 
As  indicated  in  (b)  the  structural  transformation  consists  of 
stacking-fault  formations  and  slip  (by  three  atomic  rows). 


stacking  faults  a  slip  region,  where  layers  15-18  shifted 
by  three  atomic  rows  to  an  equivalent  registry,  is  detect¬ 
ed.  As  seen  the  formation  of  the  stacking  fault  involves 
layers  12  and  13.  The  details  of  this  structural  transfor¬ 
mation  are  shown  in  Fig.  5  where  particle  positions  in 
layers  12  (solid  circles)  and  13  (open  circles)  at  three 
times  [before  (a),  during  (b),  and  after  (c)],  separated  by 
170  At,  are  shown. 

The  mechanical  response  and  energetic  characteristics 
of  the  system  can  be  investigated  best  via  layer  decompo¬ 
sition  of  the  system  properties.  In  Figs.  6(a)  aud  6(b)  wc 
depict  the  per-layer  xz  component  of  the  internal  stress 
for  the  interface  layers  (layer  5,  8,  and  9  of  the  hard  ma¬ 
terial  and  layers  10-14  of  the  soft  material).  From  these 
figures  we  observe  that  the  generation  of  the  structural 
transition  involves  a  gradual  decrease  in  the  internal  a  „ 
component  in  layers  1 1-14  of  the  soft  material  while  the 
variation  in  the  stress  in  layer  10  [the  atomic  plane  of  the 
soft  materia]  ( B )  adjacent  to  the  hard  materia]  (A),  see 
Fig.  1]  is  smaller.  As  seen  from  Fig.  6(a)  the  stress  relief 
process  and  the  associated  structural  transformation  in 
the  soft  material  are  accompanied  by  a  stress  accumula¬ 
tion  in  the  interfacial  region  of  the  hard  material  (layer 
9).  In  addition  we  observe  periodic  oscillations  in  the 
internal  stress  (most  pronounced  for  layers  11-14)  past 
the  structural  transformation  period,  as  the  system  re¬ 
laxes  in  the  new  state  after  the  structural  transformation 
events.  The  energetics  of  the  system  is  explored  via  the 
time  records  of  the  per-layer  potential  energies  (£/ )  and 
temperature  ( Tf )  shown  in  Figs.  7(a)  and  7(b),  respective¬ 
ly.  From  the  potential-energy  curves  we  find  that  the  po¬ 
tential  energy  of  particles  in  layer  1 1  (second  layer  from 
the  interface,  see  Fig.  1)  in  the  soft  material  is  initially 
lower  than  that  of  layers  12-14,  which  are  further  re¬ 
moved  from  the  interface,  since  particles  in  that  layer  are 
within  the  range  of  the  stronger  interaction  with  the  hard 
substrate  [e  ab  =s^eAAeBB '  ]•  The  potential  energy  of 

particles  in  the  interfacial  soft-material  layer  (layer  10) 
adjacent  to  the  hard  material  is  lower  by  about  —  0.4e 
than  the  value  for  layer  1 1,  and  as  seen  from  Fig.  6^a)  the 
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FIG.  5.  Details  of  the  system  trajectories  leading  to  the  for¬ 
mation  of  a  stacking  fault.  Atoms  in  layers  12  (solid  circles)  and 
13  (open  circles),  in  the  soft  material  are  shown  before  (a),  mid¬ 
way  through  the  stacking-fault  formation  (b),  and  after  the  lay¬ 
ers  have  moved  to  the  new  registry  (c).  The  atomic  trajectories 
are  for  a  simulation  at  r„  =0.81  (as  in  Figs.  3  and  4). 


extra  stabilization  “pins”  this  interfacial  layer  to  the  hard 
substrate.  The  potential  energy  of  the  topmost  layer  of 
the  hard  material  [layer  (9)]  is  lowered  further  by  about 
— 2.6e.  The  potential  energies  of  both  layers  9  and  10  in¬ 
crease  due  to  the  structural  transformation.  Since  the 
system  in  this  set  of  calculations  is  thermally  isolated,  the 
structural  change  is  accompanied  by  a  temperature  in¬ 
crease  as  seen  from  Fig.  7(b)  (where  the  curves  for  succes¬ 
sive  layers,  starting  from  10  and  up,  are  displaced  verti¬ 
cally  by  0.04e).  Note  however  that  the  final  temperature 
after  the  transformation  is  below  the  melting  temperature 
for  both  materials. 

The  variations  of  the  external  stress  components  for 
the  total  system  [calculated  from  Eq.  (13),  for  t„=0.81 
and  using  the  dynamically  determined  values  of  the  ma¬ 
trix  £f]  are  displayed  in  Fig.  8,  versus  time.  We  observe 
an  increase  in  all  components  at  the  time  of  the  structural 


FIG.  6.  Per-layer  internal  stresses  io'a )  vs  time  (in  t.u.)  for  a 
simulation  at  t  =0.81  under  adiabatic  conditions.  Layers  5,  8, 
and  9  in  the  hard  material  and  10  and  11  in  the  soft  one  are 
shown  in  (a).  The  stresses  in  layers  12-14  in  the  soft  material 
are  shown  in  (b).  The  interface  is  between  layers  9  and  10.  Note 
the  variation  in  the  internal  stresses  at  the  time  of  the  structural 
transformation.  The  internal  stresses  in  the  interior  of  the  soft 
materia)  (10-14)  decrease  with  layer  10  exhibiting  pinning  by 
the  hard  material.  The  stress  relief  in  the  soft  material  is  ac¬ 
companied  by  stress  accumulation  in  the  hard  material  (layers 
5-9). 


transformation.  We  also  note  the  axial  component,  au, 
in  the  <.  direction  (along  the  [111]  direction,  i.e.,  normal 
to  the  (111)  planes)  changes  character  from  compressive 
(due  to  the  initial  load  on  the  system)  to  tensile  (positive 
value)  past  the  structural  transformation. 

2.  Isothermal  conditions 

To  investigate  the  effect  of  the  ambient  thermal  condi¬ 
tions  we  repeated  the  simulations  described  in  (1)  but 
controlled  the  temperature  of  the  system,  at  a  reduced 
temperature  of  0. 1 1 ,  via  periodic  scaling  of  particle  veloc¬ 
ities.  This  system  also  exhibited  structural  changes  (in¬ 
terlayer  slip  and  stacking-fault  generation,  but  at  a  higher 
value  of  t,  =0.84,  as  compared  to  the  value  for  the  adia¬ 
batic  case  (r,  =0. 80).  In  addition,  the  volume  increase 
after  the  transformation  was  very  small.  Applying  higher 
values  of  r  to  the  system  it  eventually  yielded  at  a  critical 
value  rf =0.95(±0.01 ),  again  higher  than  in  the  adiabat¬ 
ic  case  (rc  =0. 86±0.02,  see  Table  I).  We  conclude  that 
the  isothermal  system  withstands  higher  values  of  the 
external  perturbation  as  compared  to  the  adiabatic  case. 
Since  the  heat  generated  during  the  transformation, 
which  in  the  adiabatic  case  can  be  utilized  to  overcome 
potential  barriers  for  structural  transformations  and 
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FIG.  7.  Per-layer  potential  energies,  E^,  in  (a),  and  tempera¬ 
tures,  T',  in  (b),  for  the  system  simulated  at  ra  *0.81  under  adi¬ 
abatic  conditions. 


yield,  is  dissipated  to  the  reservoir  under  isothermal  con¬ 
ditions,  larger  external  forces  are  required  in  the  latter 
case  in  order  to  bring  about  similar  effects. 

The  behavior  versus  time  of  the  xz  component  of  the 
internal  stress  (au )  in  layers  9- 14  at  an  applied  r=0. 87, 
is  shown  in  Fig.  9.  It  is  seen  that  while  the  stresses  in  lay¬ 
ers  10- 14  of  the  soft  material  vary  during  the  structural 
transformation  they  settle  to  values  close  to  the  initial 
ones.  Inspection  of  the  real-space  particle  trajectories  for 
the  system  reveals  that  the  system  does  evolve  through 
deformed  structures  which  contain  stacking  faults  but  it 
does  not  stabilize  in  these  configurations.  Furthermore, 
in  the  final  state  the  intralayer  registry  is  intact,  which  is 
consistent  with  the  behavior  of  the  internal  stresses 
shown  in  Fig.  9. 

3.  The  tffect  of  atomic-size  mismatch 

In  the  simulations  discussed  above  the  materials  on 
both  sides  of  the  interface  differed  by  the  potential  well- 
depth  parameter  e  [fBB  =Q.5fAA  ),  but  were  character- 


(13)]  for  simulations  at  r„*0.81  under  adiabatic  conditions  vs 
time  (in  t.u.).  Increases  are  observed  in  all  components  at  the 
time  of  structural  transformation.  Note  the  change  from 
comprehensive  (due  to  the  initial  load)  to  tensile  character  in 
the  axial  component  a"1. 


FIG.  9.  Per-layer  internal  stresses  (a1,,)  vs  time  (in  t.u.)  for  a 
simulation  under  isothermal  conditions  at  r„  =0.87.  (a)  Layers 
9-11;  (b)  layers  12-14. 
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ized  by  the  same  value  of  the  “atomic-size”  parameter, 
° aa  =a bb=(J Ab  [see  Eq.  (30)].  To  investigate  the 
dependence  of  the  system  properties  on  the  mismatch  in 
atomic  sizes  between  the  two  interfacing  materials  we 
have  performed  simulations  in  which  the  well-depth  pa¬ 
rameters  were  chosen  as  before  (see  Sec.  Ill  A)  but  the  a 
parameters  were  chosen  such  that  the  atoms  in  the  soft 
( B )  material  are  characterized  as  having  a  larger  size,  i.e., 
oBB  =  1.50^  and  a  AB  =  1 . 25cr  ^  A.  In  these  simulations 
the  hard  material  ( A )  occupied  six  layers  (numbered 
1-6)  and  the  soft  material  occupied  nine  layers  (num¬ 
bered  7-15),  with  162  atoms  per  layer  in  the  A  material 
and  72  atoms  per  layer  in  the  B  (soft)  material  (the  total 
number  of  particles  was  1620). 

A  series  of  studies,  similar  to  those  described  above, 
employing  the  {9i,r,N)  ensemble  under  thermally  adia¬ 
batic  conditions  were  performed.  From  these  simulations 
we  have  determined  that  the  critical-yield  value  of  the 
thermodynamic  tension  in  this  system  is  rexz=0.62 
±0.01  (and  the  corresponding  external  stress  aKX2  =0.51 
±0.01),  which  are  considerably  lower  than  the  corre¬ 
sponding  values  found  for  the  equal-atomic  size  systems 
(see  Table  I).  Moreover,  unlike  the  previous  cases,  where 
yield  was  preceded  by  an  external  stress  regime  in  which 
the  system  responded  inelastically  via  a  sequence  of 
structural  transformations,  for  the  present  system  (in 
which  the  atomic  sizes  across  the  interface  differ)  no  such 
behavior  is  observed.  Time  histories  of  the  layer  decom¬ 
posed  potential  energy,  Ej,  and  internal  stress,  ol„,  ob¬ 
tained  in  a  simulation  at  a  constant  thermodynamic  ten¬ 
sion  r„  =0.52  are  shown  in  Figs.  10  and  11.  First  we  ob¬ 
serve  that  at  yield  the  potential  energy  and  internal  stress 
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FIG.  10.  Per-layer  potential  energies  for  a  system  containing 
an  interface  (between  layers  6  and  7)  between  a  hard  (layers 
1-6)  and  soft  (layers  7-15)  materials  which  are  characterized  in 
addition  by  different  atomic-size  parameters  (oBB  =  \.5oAA). 
Results  are  for  simulations  employing  the  0H,t,N)  ensemble, 
under  adiabatic  conditions  at  rxt  =0.52.  Time  is  in  units  of  t.u. 
and  potential  energy  in  reduced  Lennard-Jones  units.  Note  that 
the  whole  soft  material  responds  in  unison  (layers  7-10).  The 
increase  in  potential  energy  starting  at  — 172  t.u.  corresponds  to 
yield.  No  distinct  structural  transformations  prior  to  yield  (in 
contrast  to  the  corresponding  equal-atomic-size  case  discussed 
previously)  were  detected. 
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FIG.  11.  Per-layer  internal  stresses  in  the  system  described  in 
Fig.  10.  Note  the  stress  relief,  resulting  in  yield  indicated  by  the 
vanishing  of  the  internal  stresses. 


vanish  (stress  relief).  Secondly  we  note  that  the  soft  sys¬ 
tem  (layer  numbers  larger  than  6)  does  not  exhibit  “pin¬ 
ning”  by  the  underlying  hard  substrate,  although  the 
strength  of  interaction  between  the  two  is  the  same  as  in 
the  previous  studies  where  such  pinning  was  observed  at 
the  interface  (see  above).  Thirdly  the  whole  soft  system 
responds  in  unison,  starting  at  the  interfacial  layer  and 
up  into  the  material.  These  observations  can  be  under¬ 
stood  when  considering  that  as  a  consequence  of  the 
larger  atomic  size  the  atoms  of  the  soft  material  at  the  in¬ 
terface  average  over  the  corrugation  of  the  potential  due 
to  the  substrate,  resulting  in  an  effective  potential  surface 
which  exhibits  smaller  variations  for  lateral  displace¬ 
ments  parallel  to  the  interface  plane,  and  consequently  a 
reduced  resistance  to  shear. 

C.  Applied  strain  simulations 

We  turn  next  to  another  mode  of  simulation  in  which 
instead  of  applying  a  thermodynamic  tension,  x,  the  sys¬ 
tem  is  strained  while  maintaining  a  constant  volume. 
The  preparation  of  the  system  up  to  the  application  of 
external  forces  is  the  same  as  described  in  Secs.  Ill  A  and 
III  B.  The  starting  point  for  these  simulations  is  an  equi¬ 
librium  averaged  system  of  18  layers  consisting  of  nine 
layers  each  of  hard  and  soft  materials  (see  Sec.  Ill  A),  at  a 
reduced  temperature  T  =0. 11  under  a  load  of  0.5.  (The 
equilibrium  averaged  calculational  cell  is  nearly  orthogo¬ 
nal  with  a  parallel  piped  base  with  =  1 1.04a, 
Hyy  =  6.69cr,  Hn  —  16.25o,  #^  =  3.860,  and  all  other 
components  fluctuating  about  zero.)  At  this  stage  we  ap¬ 
ply  a  constant-strain  rate  £  in  the  (112)  plane  along  the 
[110]  direction  (see  Fig.  1),  i.e., 

#o0=O  (for  a/3  for  which  a^x  and  fi^z) 

with  the  strain  rate  y=4x  10-3  t.u.-1  (i.e.,  at  this  strain 
rate  the  value  of  Hu  grows  from  zero  to  about  30%  of 
Hu  in  104  integration  time  steps). 
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As  in  the  previous  simulations  (see  Sec.  Ill  B)  we  have 
performed  studies  under  both  thermally  adiabatic  and 
isothermal  conditions.  Since  the  results  are  similar  in 
both  cases  we  focus  on  the  former  ones.  Variations  of  the 
system  temperature,  T,  and  particle  potential  energy,  Ep , 
versus  the  strain  y  are  shown  in  Fig.  12.  As  seen  these 
quantities  exhibit  a  characteristic  nonmonotonous  behav¬ 
ior  with  an  increasing  trend.  The  temperature  of  the  sys¬ 
tem  grows  in  a  stepwise  manner  and  the  potential  energy 
in  a  sawtooth  fashion  characterized  by  periods  of  increase 
in  Ep  followed  by  sharp  drops  occurring  at  the  same 
strain  values  for  which  the  step  increases  in  T  occur.  To 
assess  possible  system  size  dependencies  of  the  results  we 
display  in  Fig.  13  similar  results  for  a  system  in  which  the 
number  of  hard  material  layers  (1-9)  is  as  above  but  the 
number  of  layers  of  the  soft  material  (eBB=eAA  /2)  is 
doubled  (layers  10-27).  This  system  was  equilibrated  as 
before  and  was  simulated  at  a  constant-strain  rate  of 
2.8  X  10“ 3  t.u.-1  (which  is  equivalent  to  that  used  for  the 
smaller  system).  As  evident  the  overall  characteristics  of 
the  system  temperature  and  potential  energy  are  similar 
in  both  cases,  particularly  for  strains  up  to  and  including 
the  first  nonmonotonic  feature  (y=  0  to  ~0. 15).  Com¬ 
parison  of  the  systems’  response  at  higher  strains  indi¬ 
cates  that  the  presence  of  a  larger  region  of  soft  material 
makes  the  system  more  pliant.  However  as  we  discuss 
below  the  overall  picture  and  certain  quantitative  con¬ 
clusions  are  the  same  for  the  two  system  sizes. 

To  gain  further  detail  about  the  process  underlying  the 
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FIG.  13.  Same  as  Fig.  12,  but  for  a  larger  system  size;  nine 
layers  of  hard  material  (layers  1-9)  and  18  layers  of  soft  materi¬ 
al  (layers  10-27). 


behavior  exhibited  in  Figs.  12  and  13,  we  display  in  Figs. 
14  and  15  temperature  profiles,  T1,  decomposed  into  lay¬ 
ers  (and  offset  between  layers  by  0.03c  on  the  temperature 
axis)  versus  strain  (y)  for  the  18-  and  27-layer  systems, 
respectively.  As  seen  the  layer  temperatures  fluctuate  in 


FIG.  12.  Temperature  ( T)  and  potential  energy  ( Ep )  vs  strain 
( y )  for  a  system  composed  of  9  layers  of  hard  material  and  9 
layers  of  soft  material  (fM=f^/ 2),  simulated  under 
constant-strain  rate  (y  =  4x  10" 3  t.u.'1)  and  thermally  adiabat¬ 
ic  conditions.  The  steps  in  the  curves  are  associated  with 
structural  transformations.  All  quantities  are  given  in  reduced 
Lennard-Jones  units. 


7 

FIG.  14.  Per-layer  temperatures,  T1,  corresponding  to  the 
system  described  in  Fig.  12.  Note  the  general  heating  trend 
(thermally  adiabatic  simulations)  and  the  stages  of  sharp  tem¬ 
perature  increases,  occurring  for  strain  y  £  0.08,  in  layers 
14-18,  in  the  soft  material,  corresponding  to  structural  trans¬ 
formations. 
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FIG.  15.  Same  as  Fig.  14,  for  the  larger  system  size  (nine  lay* 
ers  of  hard  and  18  layers  of  soft  materials).  Note  sharp  temper¬ 
ature  increases  in  layers  19-27. 


response  to  the  applied  strain  with  an  overall  heating 
trend.  We  remind  the  reader  that  in  the  18-layer  system 
the  interface  between  the  hard  and  soft  material  is  be¬ 
tween  layers  9  and  10  (and  due  to  the  periodic  boundary 
conditions  also  between  1  and  18),  and  in  the  27-layer 
system  it  is  located  between  layers  9  and  10  (  and  also  be¬ 
tween  layers  1  and  27).  A  more  detailed  inspection  of 
Figs.  14  and  15  reveals  sharp  temperature  spikes  in  re¬ 
gions  inside  the  soft  materials  (layers  14-18  in  Fig.  14 
and  layers  19-27  in  Fig.  15)  occurring  at  characteristic 
values  of  the  strain  which  compare  with  the  values  of  y 
at  which  the  nonmonotonic  behavior  is  seen  in  Figs.  12 
and  13.  These  observations  indicate  that  in  response  to 
the  applied  strain  the  system  deforms  by  going  over  po¬ 
tential  barriers  between  structural  configurations.  In  fact 
examination  of  the  real-space  particle  trajectories  shows 
that  at  these  characteristic  strains  interlayer  slip  is  ini¬ 
tiated  between  layers  15  and  16  for  the  18-layer  system 
and  layers  23  and  24  for  the  larger  one  (both  inside  the 
soft  material  in  regions  removed  from  the  material  inter¬ 
face  by  about  3-5  atomic  layers). 

To  gain  further  insight  about  the  process  of  deforma¬ 
tion  we  show  in  Figs.  16  and  17  layer  profiles  for  the  per- 
layer  internal  stresses  versus  strain  in  the  18-  and  27 -layer 
systems,  respectively.  It  is  seen  that  the  stress  com¬ 
ponent  parallel  to  the  strain  direction,  a„,  exhibits  a 
sawtooth  shape  (much  like  the  potential  energies  shown 
in  Figs.  12  and  13)  with  the  first  peak  occurring  at 
Hu—0.0SHa  (i.e.,  at  a  strain  y  =0.08).  The  rises  in  the 
stress  correspond  to  stress  accumulation  and  the  subse¬ 
quent  drops  to  stress  relief  (as  corroborated  by  inspection 
of  real-space  particle  trajectories  which  reveal  that  the 
stress  relief  mechanism  involves  interlayer  slip).  The  os- 


FIG.  16.  Per-layer  internal  stress  vs  strain  ( y )  profiles  for  the 
tr'a  [in  (a)],  a'n  [in  (b)J,  and  a'^  [in  (c)J,  components  for  the  sys¬ 
tem  composed  of  nine  layers  of  hard  and  nine  layers  of  soft  ma¬ 
terial,  and  simulated  under  adiabatic  conditions.  Note  the 
monotonic  increases  in  internal  stress  in  a'a  and  the  axial  com¬ 
ponent  a'u  followed  by  sharp  drops  corresponding  to  structural 
transformations.  The  solid  dots  in  (a)  give  the  values  of  the 
internal  stress  component  o'a  obtained  via  constant  strain  simu¬ 
lations  at  the  corresponding  values  of  y.  As  seen,  for  values  of 
y  <0.063  (first  two  dots)  the  values  of  o'u  thus  obtained  coin¬ 
cide  with  those  corresponding  to  the  constant  strain  rate  simu¬ 
lations.  For  y= 0.063  ±0.002  the  internal  stress  drops  sharply, 
corresponding  to  a  structural  stress  relieving  transformation  in 
the  system.  The  constant-strain  simulations  were  used  to  deter¬ 
mine  the  critical  stresses  given  in  Table  I.  All  quantities  are  in 
reduced  Lennard-Jones  units. 


dilatory  behavior  seen  in  Figs.  16(a)  and  17(a)  after  the 
slip  is  due  to  damped  planar  oscillations  in  the  new 
structural  configuration.  The  same  general  behavior  is 
seen  for  the  axial  stress  component,  cra  (i.e.,  the  stress 
component  along  the  [111]  direction  [see  Figs.  16(b)  and 
17(b)].  Note  that  for  au  the  curves  are  bunched  in 
groups  with  the  interfacial  layers  (8,10,  and  9,11)  bound¬ 
ing  the  closely-packed  stress  curves  corresponding  to  the 
middle  region  in  the  hard  and  soft  materials  (layers  6,7 
and  12,13  for  both  system  sizes).  Finally,  the  axial  stress 
component  <7„,  in  the  direction  [llO]  of  the  applied 
strain,  shown  in  Figs.  16(c)  and  17(c),  exhibits  merely 
fluctuations  with  a  rising  trend  [also  seen  for  aa,  see 
Figs.  16(b)  and  17(b)]  which  is  a  direct  consequence  of  the 
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FIG.  17.  Same  as  Fig.  16,  but  for  the  larger  system,  com¬ 
posed  of  nine  layers  of  hard  material  and  18  layers  (10-28)  of 
soft  material. 


heating  of  the  system.  Note  in  addition  the  transition 
from  compressive  (due  to  the  initial  load)  to  tensile  be¬ 
havior  as  the  strain  increases.  The  difference  between  the 
character  of  variations  in  (similarly  a„ )  and  that  of 
au  and  an  is  consistent  with  our  observation  that  (he 
stress  relief  mechanism  is  via  processes  in  which  (111) 
planes  slip  past  one  another  with  no  significant  intralayer 
distortions. 

As  already  mentioned  in  Sec.  IIIB  the  rate  of  applica¬ 
tion  of  the  perturbation  [applied  shear  stress  (or  thermo¬ 
dynamic  tension)  as  discussed  in  Sec.  IIIB  and  applied 
strain  in  this  section]  in  MD  simulations  is  faster  than 
laboratory  rates.  In  addition  these  rates  are  also  too  fast 
to  allow  an  adiabatic  (in  the  sense  of  dynamical  relaxa¬ 
tion  on  the  atomic  scale)  response  of  the  system.  With  re¬ 
gard  to  the  first  issue  we  refer  to  the  studies  of  Moran, 
Ladd,  and  Hoover2*  who  suggested  on  the  basis  of  com¬ 
parison  between  NEMD  simulations  and  experiments 
that  large-deformation  and  plastic-flow  phenomena  in 
close-packed  metals  obey  a  scaling  relation  which  allows 
comparison  between  MD  calculations  and  experimental 
data  although  the  two  differ  by  orders  of  magnitude  in 
the  rates  of  the  applied  shear.  Moreover  these  authors 
suggest  that  these  observations  imply  that  plastic 
response  and  flow  of  these  metals  can  be  described  by  a 
single  physical  mechanism  over  a  range  of  strain  rates 
from  10  kHz  to  1  THz  (our  values  of  y  using  well  depth, 
e,  and  a  parameter  values  in  the  6-12  LJ  potential  to  ap¬ 


proximate  typical  metals,  are  of  the  order  of  10  GHz). 

To  address  the  second  issue  and  determine  accurately 
temporally  adiabatic  values  for  the  critical  strains  (i.e., 
those  at  which  inelastic  structural  deformations  occur)  in 
our  system,  we  have  performed  for  both  system  sizes 
simulations  at  constant  strains  approaching  the  value  at 
which  slip  occurred  [see  first  peak  in  Figs.  16(a)  and  17(a) 
at  y — 0.08].  Thus  in  these  series  of  simulations  the  sys¬ 
tem  was  allowed  to  fully  relax  under  prescribed  strains. 
The  values  for  the  internal  stresses,  cr^,  obtained  from 
these  simulations  are  represented  by  solid  dots  in  Figs. 
16(a)  and  17(a).  First  we  note  that  for  strains  below  the 
critical  strain  value  the  developed  internal  stresses  in 
these  simulations  are  the  same  as  those  obtained  in  the 
constant-strain-rate  studies,  which  indicates  that  for  this 
regime  of  strain  values  the  system  responds  adiabatically 
even  when  the  strain  is  applied  at  a  rather  high  rate. 
Secondly,  the  values  of  the  critical  strain,  ys  obtained 
in  these  simulations  (see  Table  I,  =7/^  /Ha 
=0.063  ±0.002  for  the  18-layer  system  and  0.076±0.002 
for  the  27-layer  system)  are  below  that  given  by  the 
constant-strain-rate  simulations.  Thirdly,  the  corre¬ 
sponding  values  of  the  internal  stresses  o^a  (0.975  and 
1.14  for  the  small  and  large  system,  respectively,  see 
Table  I)  are  larger  than  the  one  obtained  using  the  adia¬ 
batic  simulations  (Sec.  IIIB),  and  the  one  for 

the  27-layer  system  is  somewhat  higher  than  that  of  the 
smaller  system.  Note  however,  that  the  critical  strain- 
to-stress  ratio,  y,  /a ,  is  nearly  the  same  for  both  systems 
(0.0634  and  0.0652,  for  the  small  and  large  systems,  re¬ 
spectively).  We  conclude  therefore  that  while  the  larger 
system  exhibits  a  somewhat  larger  tolerance  to  strain 
(i.e.,  is  more  pliant)  the  elastic  properties  of  our  systems 
are  independent  of  thickness  at  or  above  nine  layers  of 
soft  material.  Finally,  we  note  that  for  strain  values 
above  y„  the  systems  exhibit  slip  accompanied  by  sharp 
drops  in  the  value  of  the  internal  stress  a„  (and  similarly 
for  o a )  indicating  stress  relief.  We  note  that  removing 
the  strain  for  strain  values  below  the  critical  value  ys,  re¬ 
sults  in  a  return  of  the  system  to  the  unstrained 
configuration,  indicating  that  in  this  strain  regime  the 
system  deforms  elastically,  while  after  removing  the 
strain  past  the  structural  transformation  (i.e.,  for  values 
larger  than  y, )  the  system  remains  in  the  new  structural 
configurations. 


IV.  SUMMARY 

In  this  investigation  we  have  studied  the  formulations 
of  molecular-dynamics  techniques  for  simulations  of  ma¬ 
terial  systems  evolving  under  applied  finite  external  per¬ 
turbations.  Methodologies  of  simulations  of  phenomena 
involving  finite  deformations  (plastic  deformations, 
structural  transformations,  yield,  and  flow)  have  been  dis¬ 
cussed  (Sec.  II)  and  demonstrated  (Sec.  III).  These  simu¬ 
lations  yield  valuable  information  about  the  atomic  scale 
mechanisms  of  materials’  dynamical  response  to  mechan¬ 
ical  perturbations. 

In  these  studies  we  focused  on  interfacial  systems,  com¬ 
posed  of  interfacing  crystalline  solids  characterized  by 
differing  interatomic  interaction  strengths  and  in  some 
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cases  also  by  differing  atomic  sizes.  Our  main  results 
may  be  summarized  as  follows. 

(i)  For  interfacial  systems  which  are  characterized  by 
differing  interatomic  interaction  strengths  (i.e.,  the  inter¬ 
face  is  between  a  hard  and  soft  material),  the  system 
responds  to  an  applied  nonisotropic  perturbation  (applied 
shear  or  strain)  first  elastically  and  then  via  stress  relief 
mechanisms  which  involve  structural  transformations 
(stacking-fault  formation  and  interlayer  slip).  For  larger 
values  of  the  external  forces,  eventual  yield  occurs. 

(ii)  Critical  values  of  the  external  perturbation  required 
in  order  to  bring  about  inelastic  response  (structural 
transformations  and  yield)  have  been  determined  (see 
Table  I).  Our  simulations  demonstrate  that  these  critical 
values  are  smaller  for  a  system  under  thermally  adiabatic 
conditions  than  for  an  isothermal  environment.  The  ori¬ 
gin  of  the  difference  lies  in  the  dissipation  of  the  generat¬ 
ed  thermal  energy  under  isothermal  conditions,  which 
necessitates  larger  values  of  the  external  perturbations  in 
order  to  overcome  potential  barriers  for  structural 
modifications  and  eventual  yield. 

(iii)  The  cohesive  interatomic  interactions  at  the  inter¬ 
face  between  a  hard  substrate  and  a  soft  material  result  in 
"pinning”  of  the  soft  material  at  the  interface  (1-3  atom¬ 
ic  layers).  As  a  result  the  response  of  the  system  to  the 
external  perturbation  (stress  relief  via  plastic  as  well  as 
structural  transformations)  occurs  in  a  “shear  band”  con¬ 
sisting  of  a  few  atomic  layers  inside  the  soft  material, 
which  are  located  at  about  1-3  layers  away  from  the 
original  (unstressed)  interface.  The  stress  relief  in  the  soft 
material  is  accompanied  by  stress  accumulation  in  the 


hard  substrate. 

(iv)  MD  simulations  for  interfacing  hard  and  soft  ma¬ 
terials,  which  in  addition  are  characterized  by  differing 
atomic  sizes,  reveal  the  important  role  played  by  atomic- 
size  mismatch  in  determining  the  atomic-scale  mecha¬ 
nism  of  response.  For  such  system  it  was  found  that  no 
adhesive  pinning  occurs  at  the  interface,  and  that  the  soft 
(and  larger  atomic  size)  material  responds  as  a  whole  with 
no  distinct  structural  transformations  preceding  the  yield 
point.  The  critical-yield  stress  value  for  this  system  is 
significantly  lower  than  that  found  for  the  corresponding 
equal-atomic-size  system. 

(v)  Comparison  of  our  results  in  this  study  for  the  [111] 
interface  with  our  previous  investigations  of  the  [100]  in¬ 
terface11  demonstrates  the  dependence  of  the  critical 
values  of  the  shear  stresses  on  the  crystallographical 
orientation  of  the  interface,  as  well  as  of  certain  details  of 
the  stress  relief  mechanisms. 
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35If  the  interaction  parameters  are  chosen  such  that  they  corre¬ 
spond  to  the  cohesive  energy  and  lattice  constant  of  nickel 
(e=3.54x  10-11  erg,  <r=2.49  A,  and  atomic  mass 
m  =9.75x  10~23  g)  a  reduced  temperature  T  =0. 1 1  corre¬ 
sponds  to  300  K,  the  reduced  melting  temperature  Tm  =0.7 
corresponds  to  2000  K,  the  reduced  unit  of  stress  or  load  to 
2.4X  107  g/cm2  (or  2.4x  10'°  dyn/ctn2),  and  the  reduced  time 
unit  (tu.)  corresponds  to  4. 1 X 10“ 13  sec. 
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Basic  understanding  of  the  structure  and  dynamics  of  materials  and  their  response  to  external 
perturbations  requires  knowledge  on  the  microscopic  level,  of  the  underlying  energetics  and 
atomic  dynamics,  whose  consequences  we  observe  and  measure.  Coupled  with  the  above  is  the 
everlasting  quest  to  observe  and  understand  natural  phenomena  on  refined  microscopic  scales, 
which  provides  the  impetus  for  the  development  of  experimental  and  theoretical  techniques  for 
the  interrogation  of  materials  with  refined  spatial  and  temporal  resolution.  In  this  paper  we 
review  molecular  dynamics  simulations  for  studies  of  the  energetics  and  dynamical  response  of 
materials  to  external  mechanical  perturbations.  Applications  to  investigations  of  solid  and  liquid 
interfacial  systems  under  stress  and  to  studies  of  the  consequences  of  tip-substrate  interactions  in 
atomic  force  microscopy  are  demonstrated. 


I.  INTRODUCTION 

Basic  understanding  of  the  structure  and  dynamics  of  mate¬ 
rials  and  their  properties  often  requires  knowledge  on  a  mi¬ 
croscopic  level  of  the  underlying  energetics  and  interaction 
mechanisms  whose  consequences  we  observe  and  measure. 
The  everlasting  quest  to  observe  and  understand  nature  on 
microscopic  scales1  is  a  dominant  trend  in  science,  leading  to 
the  development  of  conceptual,  technological,  and  experi¬ 
mental  devices  which  allow  interrogation  of  nature  with  re¬ 
fined  spatial  and  temporal  resolution.  Elucidation  of  the 
structure,  dynamics,  mechanical  properties,  and  response  of 
material  systems  to  external  perturbations  on  the  atomic  lev¬ 
el,  are  key  issues  in  developing  a  fundamental  understanding 
of  varied  systems  and  phenomena  of  coupled  basic  and  tech¬ 
nological  interest,  such  as  surface  and  interface  systems,  sur¬ 
face  reactions  and  catalysis,  electrochemical  interfaces, 
understanding  of  structure-function  relationships  in  bio- 
molecular  systems,  microelectronics  materials,  tribology, 
lubrication,  wear,  material  fatigue  and  yield,  crack  propaga¬ 
tion,  stress  induced  phase  and  structural  transformations, 
and  hydrodynamics!  phenomena  in  confined  fluid  systems, 
to  name  just  a  few.  Although  most  of  the  above  listed  phe¬ 
nomena  represent  everyday  experiences  and  have  been  ob¬ 
served  and  studied  for  a  long  time,  detailed  microscopic  the¬ 
ories  of  them  (with  few  exceptions)  are  lacking. 
Nevertheless  large  bodies  of  empirical  data  and  in  some 
cases  phenomenological  model  descriptions  have  been  de¬ 
veloped. 

On  the  experimental  front  the  developments  of  the  surface 
force  apparatus  (SFA)U  scanning  tunneling  microscopy 
(STM)4  and  of  the  related  atomic  force  microscopy 
(AFM)5  broaden  our  perspectives  and  abilities  to  probe  the 
morphological  and  electronic  structure  and  the  nature  of 
interatomic  forces  in  materials,  open  new  avenues6  for  mi¬ 
croscopic  investigations  and  manipulations  of  technological 
systems  and  phenomena  such  as  tribology,3,7  lithography 
and  in  biochemical  applications,6  and  present  exciting  possi¬ 
bilities  for  the  development  of  miniaturized  electronic,  mag¬ 
netic  and  optical  devices. 


In  the  latter  two  techniques  a  sharp  tip  is  brought  close  to 
the  surface  and  either  the  tunneling  current  (STM)  or  de¬ 
flection  of  a  flexible  cantilever  holding  the  tip  (AFM)  are 
monitored,  either  by  tunneling  or  by  optical  interferometry 
as  the  surface  is  scanned.  Since  the  proximity  of  the  tip  to  the 
surface  may  induce  structural  modifications  in  both  the  sub¬ 
strate  and  the  tip,  the  nature  of  the  tip-substrate  interactions 
and  their  consequences  are  key  issues  in  the  development  of 
these  spectroscopies  (as  is  the  issue  of  probe-system  interac¬ 
tion  in  all  methods  of  microscopic  resolution).  Of  particular 
interest  are  questions  related  to  the  dynamical  response  of 
the  substrate  (and  tip)  which  may  result  in  temporary  or 
permanent  modifications  of  the  local  properties,  and  those 
related  to  the  effect  of  system  conditions  (such  as  tempera¬ 
ture)  which  may  influence  the  response  as  well  as  affect  the 
spatial  resolution.  It  has  been  indeed  suggested1*  that  STM 
images  can  be  dominated  by  elastic  deformations  induced  by 
the  tip-substrate  interatomic  forces  and  a  transition  from  a 
tunnel  regime  to  point  contact  has  been  observed.’ 

Recently,  the  AFM  has  been  applied  to  study  the  topo¬ 
graphy  and  frictional  forces  at  the  atomic  scale.7  Further¬ 
more,  recent  investigations3  of  the  dynamical  properties  of 
molecularly  thin  liquid  films  confined  between  two  surfaces, 
using  an  extension  of  the  SFA,  have  demonstrated  the  direct 
consequences  of  microscopic  molecular  structure  and  inter¬ 
actions  on  the  structural  and  dynamical  properties  of  such 
sytems  both  under  static  and  steady-state  shear-flow  condi¬ 
tions. 

On  the  theoretical  front  computer  molecular  dynamics 
(MD)  simulations’*3"1’,  which  are  in  a  sense  computer  ex¬ 
periments  where  the  evolution  of  a  physical  system  is  simu¬ 
lated,  with  refined  temporal  and  spatial  resolution,  via  a  di¬ 
rect  numerical  solution  of  the  model  equations  of  motion, 
open  new  avenues  in  investigations  of  the  microscopic  ori¬ 
gins  of  material  phenomena.  These  methods  alleviate  certain 
of  the  major  difficulties  which  hamper  other  theoretical  ap¬ 
proaches,  particularly  for  complex  systems  such  as  those 
characterized  by  a  large  number  of  degrees  of  freedom,  lack 
of  symmetry,  nonlinearities  and  complicated  interactions.  In 
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addition  to  comparisons  with  experimental  data  which  pro¬ 
vide  guidance  in  the  interpretation  of  data  and  the  impetus 
for  new  measurements,  computer  simulations  can  be  used  as 
a  source  of  physical  information  which  may  not  be  accessible 
to  laboratory  experiments,  and  in  some  instances  the  com¬ 
puter  experiment  itself  serves  as  a  testing  ground  for  theo¬ 
retical  concepts. 

In  this  paper  we  review  recent  molecular  dynamics  studies 
of  systems  and  phenomena  related  to  the  rapidly  evolving 
field  of  nanoscience.  We  address  and  demonstrate  investiga¬ 
tions  of  three  issues:  (i)  the  consequences  of  dynamical  tip- 
substrate  interactions,  (ii)  the  nature  of  bounded  thin  fluid 
films  in  equilibrium  and  in  steady-state  flow  under  shear, 
and  (iii)  the  microscopic  energetics,  dynamics  and  response 
of  interfacial  systems  to  external  perturbations. 

II.  CONSEQUENCES  OF  TIP-SUBSTRATE 
INTERACTIONS 

To  investigate  the  consequences  of  the  dynamical  interac¬ 
tions  between  the  tip  and  the  substrate  in  AFM  we  have 
performed  MD  simulations,  employing  the  Stillinger- We¬ 
ber  (SW)  potentials  for  silicon20  which  include  two-  and 
three-body  interactions  and  which  have  been  used  recently 
in  a  number  of  investigations  of  bulk,  surface  and  interface 
properties  of  this  material. 17-24  Since  both  the  substrate  and 
tip  consist  of  the  same  material,  these  simulations  corre¬ 
spond  to  the  case  of  a  reactive  tip-substrate  system.  Our 
simulations, 17-19,23  in  both  the  constant-tip-height  and  con¬ 
stant-force  scan  modes,  reveal  that  the  local  structure  of  the 
surface  can  be  stressed  and  modified  as  a  consequence  of  the 
tip-substrate  dynamical  interaction,  even  at  tip-substrate 
separations  which  correspond  to  weak  interaction.  For  large 
separations  these  perturbations  anneal  upon  advancement  of 
the  tip  while  permanent  damage  can  occur  for  smaller  sepa¬ 
rations.  For  die  material  that  we  simulated  (Si),  we  do  not 
find  long-range  elastic  deformations,  which  may  occur  in 
other  circumstances*  depending  upon  the  elastic  properties 
of  the  material  and  the  nature  of  interactions.  Furthermore, 
we  find  that  the  characteristics  of  the  data  depend  upon  the 
geometry  of  the  scan,  the  degree  of  perfection  of  the  sub¬ 
strate  and  the  temperature.  We  identify  various  dynamical 
events  including  stick-slip  phenomena,  which  could  be  ex¬ 
perimentally  resolved,  using  current  estimates,3,6  and  which 
would  influence  the  analysis  of  data,  as  well  as  pointing  to 
ways  in  which  the  tip  could  be  used  for  atomic  scale  manipu¬ 
lations  of  the  material. 

In  our  simulations  the  system  consists  of  NL  layers  of 
dynamic  Si  particles  with  nL  atoms  per  layer,  exposing  the 
(111)  surface,  and  interacting  with  two  layers  of  a  static  Si 
substrate  of  the  same  structure.  In  simulations  involving  a 
small  tip  Nl  =  4  and  nL  =  49,  while  in  simulations  involv¬ 
ing  larger  tips  NL  =  6  and  nL  =  100.  The  two-dimensional 
calculational  cell,  defined  by  the  ( llO)  and  ( lOl)  vectors,  is 
periodically  repeated  parallel  to  the  ( 1 1 1 )  plane.  Two  types 
of  tips  were  employed  in  our  simulations.  First,  we  have  used 
a  sharp  tip  simulated  by  four  Si  atoms  in  an  initial  tetrahe¬ 
dral  configuration,  “mounted”  on  and  interacting  with  two 
layers  of  silicon  atoms  which  serve  as  a  holder.  In  addition 
we  have  employed  larger  tips  consisting  of  102  atoms,  which 


time  (tu) 

Fig.  1.  Tot*]  F,  force  (in  units  of  1.6572xl0-,  N)  on  the  small  rigid 
silicon  tip  as  it  is  lowered  to  a  distance  of  1.227  A,  equilibrated  and  scans 
along  the  surface.  Time  is  in  units  of  tu  =  7.6634  x  10” 14  s.  Note  the  sharp 
variation  in  the  force  as  the  tip  is  lowered  and  the  periodic  oscillations 
reflecting  the  atomic  structure  of  the  surface  during  the  scan  (ti.  60) . 


were  either  ordered  [exposing  a  16  atoms  (111)  planar  fac¬ 
et]  or  disordered.  The  equations  of  motion,  governed  by  the 
SW  potential,20  are  integrated  using  a  fifth-order  predictor- 
corrector  algorithm  with  a  time  step  A/  =  1.15x  10-3  ps, 
[or  0.015  time  units  (tu)  where  tu  =  7.6634  x  10- 14  s]. 
Throughout  we  use  e  =  50  kcal/mol  as  the  unit  of  energy, 
a  —  2.0951  A  as  the  unit  of  length,  and  e/ 
a  —  1.65728  X  10~9  N  as  the  unit  of  force.  The  kinetic  tem¬ 
perature  is  controlled  to  room  temperature  via  scaling  of 
particle  velocities  in  the  bottom  layer  of  the  dynamic  sub¬ 
strate.  Both  constant-tip-height  and  constant-force  scan 
modes  were  simulated. 

First  we  investigate  the  process  of  lowering  a  tip  towards 
the  substrate  which  is  then  followed  by  either  raising  of  the 
tip  or  scanning  laterally  across  the  surface  in  either  a  con¬ 
stant  height  or  constant-force  mode.  The  total  force  in  the 
normal  (z)  direction  on  a  rigid  four-atom  tip  (arranged  as  a 
tetrahedron)  is  shown  in  Fig.  1.  Starting  with  the  system  at 
equilibrium  and  outside  the  range  of  interaction  with  the 
substrate  [see  snapshot  of  the  system  initial  configuration  in 
Fig.  2(a)]  the  tip  is  lowered  slowly  toward  the  surface.  Ini¬ 
tially  the  tip  experiences  an  attractive  force  (F<0)  accom¬ 
panied  by  a  lifting  of  the  surface  atom  right  below  the  lower¬ 
most  tip  atom,  as  seen  in  the  snapshot  given  in  Fig.  2(b).  As 
the  tip  approaches  the  surface  the  attraction  diminishes  cul¬ 
minating  in  a  repulsion  ( F ,  >0  at  —30  tu  in  Fig.  1 )  which 
turns  into  an  attraction  in  a  sudden  manner  (see  Fig.  1  at 
5  40  tu ) .  Inspection  of  the  real  space  trajectory  of  the  system 
reveals  that  these  events  are  associated  with  a  formation  of 
an  interstitial  defect  where  the  surface  atom  is  pushed  into 
the  interstitial  position  in  the  second  layer.  Further  relaxa¬ 
tion  of  the  system  at  tip-to-substrate  separation  of  1.227  A 
(corresponding  to  a  repulsive  region  of  the  interatomic  po¬ 
tential)  results  in  a  net  attractive  force  caused  by  the  tip- 
induced  local  reconstruction  [see  Fig.  2(c)].  At  this  point 
the  latteral  scan  begins  (in  the  [110]  direction)  at  a  rate 
which  assures  tip-substrate  relaxed  configurations  through¬ 
out  the  scan.  As  the  tip  scans  (at  constant  height)  over  the 
surface  the  recorded  normal  force  on  it  oscillates  with  the 
periodicity  of  the  substrate,  being  repulsive  (F,  >0)  when 
the  tip  traverses  the  region  between  atoms  and  attractive 
over  the  substrate  atoms  ( Fig.  1 ,  t  £  60  tu ) .  Snapshots  of  the 
system  at  different  tip  locations  (starting  from  the  tip-sub¬ 
strate  equilibrated  configuration  shown  in  Fig.  2(d)  and 
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Fig.  2.  Snapshots  of  small  static  tip-substrate  atomic  configurations.  Tip 
atoms  in  bright  white  and  topmost  substrate  layer  as  dark  spheres,  (a)  Tip 
outside  range  of  interaction;  (b)  tip  lowered  towards  the  surface,  “lifting" 
the  atom  underneath  it;  (c)  tip  equilibrated  with  the  surface  at  the  lowered 
position  ( 1.227  A  from  the  nominal  position  of  the  surface  plane),  demon¬ 
strating  the  formation  of  an  interstitial  defect  under  the  tip  ( note  location  of 
the  dark  ball  directly  under  the  tip);  (d)  initial  cofiguration  before  scan 
(translation  of  the  tip  to  the  left)  begins;  (e)  configuration  when  the  tip 
scans  between  substrate  atoms;  ( f)  configuration  when  the  tip  translated  by 
one  unit  cell  distance,  inducing^  new  interstitial  surface  defect.  The  pre¬ 
vious  tip-induced  interstitial  [  see  ( d )  ]  returns  to  his  original  location  in  the 
first  atomic  layer. 


moving  the  tip  laterally  to  the  left)  show  that  when  in  the 
region  between  atoms  the  tip-induced  interstitial  returns  to 
his  original  location  in  the  first  atomic  layer  [see  Fig.  2(e)  ] 
and  a  new  interstitial  defect  is  generated  when  the  tip  moves 
to  a  location  over  the  next  atom  in  the  surface  row  [Fig. 
2(01- 

The  consequences  of  the  interaction  between  a  large  dy¬ 
namic  tip  consisting  of  102  atoms,  arranged  in  four  layers 
and  exposing  a  16-atom  (111)  facet,  are  exhibited  in  Fig.  3 
where  snapshots  of  the  system  configurations  initially  [Fig. 
3(a)],  at  the  lowered  tip  configuration  [Fig.  3(b)],  and 
after  raising  the  tip  [Fig.  3(c)  ],  are  shown.  As  seen,  for  this 
reactive  system,  lifting  of  the  tip  results  in  the  formation  of  a 
“connective  neck”  consisting  mostly  of  tip  atoms.  The  varia¬ 
tions  in  the  normal  component  of  the  force  on  the  tip  atoms 
Ft,  potential  energy  of  the  tip  atom  EP,  and  their  kinetic 
temperature  T,  are  shown  in  Fig.  4(b)-4(d),  respectively,  as 
well  as  the  temperature  of  the  dynamic  substrate  and  height 
of  the  tip  assembly  from  the  surface  in  Figs.  (4e)  and  4(a), 
respectively.  The  snapshots  in  Fig.  3(a)-3(c)  correspond  to 
times  t  =  0,  50,  and  120  in  Fig.  4.  As  seen  from  Fig.  4(b)  at 
the  equilibrium  configuration,  with  the  tip  resting  in  contact 
with  the  surface,  the  total  normal  component  of  the  force  on 
the  tip  atoms  is  attractive  (F,  <0  in  the  time  interval, 
20  <  /  <  60  tu ) .  Subsequent  raising  of  the  tip  ( /  >  60)  results 
in  an  increased  attraction  between  the  tip  and  the  substrate 
atoms  and  an  increase  in  the  potential  energy  [see  Fig.  4(c)  ] 
of  the  tip  atoms  as  the  initial  crystalline  structure  of  the  tip  is 
modified  due  to  strained  bonding  to  the  substrate.  Further 
raising  of  the  tip  ( t  —  1 00 )  results  in  disruption  of  some  of  the 
bonds  between  the  tip  and  substrate  atoms,  resulting  in  a 
sudden  decrease  in  the  force  between  them  [see  Fig.  4(b)  ], 


Fig.  3.  Atomic  configurations  as  a  large  faceted  dynamical  tip  initially 
equilibrated  outside  the  range  of  interaction  (a),  is  lowered  to  a  distance  of 
2.91  A  from  the  surface  (b),  and  following  equilibration  is  raised  (c). 


time  (tu) 


Fig.  4.  Variation  vs  time  (in 
tu  =  7.6634  X  10* 14  s)  of:  (a)  Center  of 
mass  height  of  the  tip-holder  assembly, 
Z;  (b)  normal  force  on  the  tip,  F,  (in 
units  of  e/<7  m  1.6572X  10_*N;  (c)  po¬ 
tential  energy  of  the  tip  atoms,  £,;  (d) 
temperature  of  the  tip  atoms,  (e)  kinetic 
temperature  of  the  dynamic  substrate. 
Simulations  are  for  lowering  and  raising 
of  a  large  (102  atom)  faceted  ordered 
dynamical  tip.  See  corresponding  atom¬ 
ic  configurations  in  Fig.  3. 
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an  increase  in  the  potential  energy  and  a  sudden  rise  in  the  tip 
and  local  substrate  temperature  [see  Figs.  4(d)  and  4(e) ], 
yielding  eventually  the  atomic  configuration  shown  in  Fig. 
3(c). 

In  constant-force  scanning  simulations,  in  addition  to  the 
particle  equations  of  motion,  the  center  of  mass  of  the 
tip-holder  jssembly,  Z,  is  required  to  obey  MZ 
=  (F(t)  —  Fm  )*Z  —  yZ  where  F  is  the  total  force  exerted 
by  the  tip-atoms  on  the  static  holder  at  time  t,  which  corre¬ 
sponds  to  the  force  acting  oja  the  tip  atoms  due  to  their  inter¬ 
action  with  the  substrate,  Fat  is  the  desired  (prescribed) 
force  for  a  given  scan,  y  is  a  damping  factor,  and  M  is  the 
mass  of  the  holder.  In  these  simulations  the  system  is 
brought  to  equilibrium  for  a  prescribed  value  of  F.„, ,  and 
the  scan  proceeds  as  described  above  while  the  height  of  the 
tip-holder  assembly  adjusts  dynamically  according  to  the 
above  feedback  mechanism. 

In  Figs.  5(a)-5(d)  and  lug.  6  we  show  the  results  for  a 
constant-force  scan,  for  F.„,  —  —  13  e/a  (corresponding 
to  —  2.15X10-8  N,  i.e.,  negative  load).  In  these  simula¬ 
tions  the  substrate  consists  of  six  dynamic  layers  with  100 
atoms  per  layer.  Side  views  of  the  system  trajectories  at  the 
beginning  and  end  stages  of  the  scan  are  shown  in  Figs.  5(a) 
and  5(b),  and  5(c),  respectively.  As  seen,  the  tip-substrate 
interactions  induce  local  modifications  of  the  substrate  and 
tip  structure,  which  are  transient  (compare  the  surface 
structure  under  the  tip  at  the  beginning  of  the  scan  [Fig. 
5(a)  ],  exhibiting  outward  atomic  displacements  of  the  top- 
layer  atoms,  to  that  at  the  end  of  the  scan  [Fig.  5(c),  where 
that  region  relaxed  to  the  unperturbed  configuration] .  The 
recorded  force  on  the  tip  holder  along  the  scan  direction  ( x ) 


<1S1>  o  t  x(«)  * 


Fig.  5.  (*)-(c)  Particle  trajectories  in  a  constant-force  simulation, 
m  —  13.0  (i.e.,  —  2.15x10“*  N),  viewed  along  the  (10T)  direction  just 
before  (a)  and  after  (b)  a  stick-slip  event  and  towards  the  end  of  the  scan 
(c),  fora  large,  initially  ordered,  dynamic  tip.  (d)  The  recorded  F,,  exhibit¬ 
ing  stick-slip  behavior,  (e)  The  FM  force  in  a  constant-force  scan  (F„„ 
•  1.0)  employing  a  glassy  static  tip,  exhibiting  the  periodicity  of  the  sub¬ 
strate.  Shown  in  the  inset  are  the  real-space  trajectories  towards  the  end  of 
the  scan,  demonstrating  the  tip-induced  substrate  local  modifications. 


Fig.  6.  Constant-force  scan  sim¬ 
ulation  at  FM1I  =  —  13  (i.e., 
negative  load  of  —  2.15x  10“K 
N),  employing  a  large  (102 
atoms)  faceted,  ordered  dynami¬ 
cal  tip.  (a)-(c)  Variations  of  the 
force  F,  (along  the  scan  direc¬ 
tion)  and  of  the  force  compo¬ 
nents  normal  to  it,  Fy  and  F, 
Forces  in  the  units  of 
1.6572x10“''  N.  Note  that  the 
recorded  force  in  the  z  direction 
fluctuates  around  the  prescribed 
value;  (d)  center-of-mass  height 
of  the  tip-holder  assembly;  (e) 
potential  energy  Ep  of  the  tip 
atoms:  (f)  kinetic  temperature  of 
the  tip  atoms.  Note  the  discon¬ 
tinuous  variation  and  asymme¬ 
try  in  F, ,  signifying  stick-slip  be¬ 
havior,  and  the  accompanying 
variations  in  E,  and  T. 


is  shown  in  Fig.  5(d)  and  in  Fig.  6(a),  exhibiting  a  periodic 
modulation,  portraying  the  periodicity  of  the  substrate.  At 
the  same  time  the  normal  force  F,  fluctuates  around  the 
prescribed  value  [Fig.  6(c)]  and  no  significant  variations 
are  observed  in  the  force  component  normal  to  the  scan  di¬ 
rection  [F,  in  Fig.  6(b)]. 

Most  significant  is  the  stick-slip  behavior  signified  by  the 
asymmetry  in  Fx  [observed  also  in  the  real-space  atomic 
trajectories  in  Figs.  5(a)  and  5(b)].  Here,  the  tip  atoms 
closest  to  the  substrate  attempt  to  remain  in  a  favorable 
bonding  environment  as  the  tip-holder  assembly  proceeds  to 
scan.  When  the  forces  on  these  atoms  due  to  the  other  tip 
atoms  exceed  the  forces  from  the  substrate,  they  move  rapid¬ 
ly  by  breaking  their  current  bonds  to  the  surface  and  forming 
new  bonds  in  a  region  translated  by  one  unit  cell  along  the 
scan  direction.  The  detailed  energetics  of  the  atomic-scale 
stick-slip  phenomenon  can  be  elucidated  from  the  variations 
in  the  potential  and  kinetic  energies  of  the  tip  atoms  along 
the  scan,  shown  in  Fig.  6(e)  and  6(0,  respectively.  As  seen, 
during  the  stick  stage,  the  potential  energy  of  the  strained 
bonds  between  the  tip  and  substrate  atoms  increases.  The 
slip  stage  is  signified  by  a  discontinuity  in  the  force  along  the 
scan  direction  [ Fx ,  in  Fig.  6(a)  ],  and  by  a  sharp  decrease  in 
the  potential  energy  [Fig.  6(e)  ],  which  is  accompanied  by  a 
sudden  increase  in  the  kinetic  temperature  of  the  tip  atoms 
[Fig.  6(0  ]  as  a  result  of  the  disruption  of  the  bonds  to  the 
substrate  and  rapid  motion  of  the  tip  atoms  to  the  new  equi¬ 
librium  positions.  We  note  that  the  excess  kinetic  energy 
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(local  heating)  acquired  by  the  tip  during  the  rapid  slip, 
dissipates  effectively  during  the  subsequent  stick  stage,  via 
the  tip  to  substrate  interaction  [  see  the  gradual  decrease  in  T 
in  Fig.  6(f),  following  the  sudden  increase]. 

We  note  that  our  constant-force  simulation  method  corre¬ 
sponds  to  the  experiments  in  Ref.  7  in  the  limit  of  a  stiff  wire 
(lever)  and  thus  the  stick-slip  phenomena  which  we  ob¬ 
served  are  a  direct  consequence  of  the  interplay  between  the 
surface  forces  between  the  tip  and  substrate  atoms  and  the 
interatomic  interactions  in  the  tip.  The  Fx  force  which  we 
record  corresponds  to  the  frictional  force.  From  the  extrema 
in  Fx  [Figs.  5(d)  or  6(a)]  and  the  load  (Fwx )  used  we 
obtain  a  coefficient  of  friction  fi  =  |FX  \/\F..„  j  =  0.77,  in 
the  range  of  typical  values  obtained  from  tribological  mea¬ 
surements  in  vacuum,  although  we  should  caution  against 
taking  this  comparison  rigorously. 

Results  for  a  constant-force  scan  at  a  positive  load  (F,„. 
=  0. 1,  i.e.,  1.66X  10“ 10  N),  employing  the  large  faceted  tip, 
are  shown  in  Fig.  7.  As  seen  from  Fig.  7(d),  the  center-of- 
mass  height  of  the  tip  and  holder  assembly  from  the  surface, 
Z,  exhibits  an  almost  monotonic  decrease  during  the  scan,  in 
order  to  keep  the  force  on  the  tip  atoms  around  the  pre¬ 
scribed  value  of  0.1  [see  Fig.  7(e)  ].  At  the  same  time  the 
potential  energy  of  the  tip  increases.  This  curious  behavior 
corresponds  to  a  “smearing”  of  the  tip  as  revealed  from  the 
real-space  trajectories  shown  in  Figs.  7(a)-7(c).  Compari- 


Fig.  7.  Constant  scan  similation  at  Fx„,  »  0.1  (i.e.,  1.6572X  10_l°  N), 
employing  a  large  (102  atoms),  initially,  ordered  dynamic  tip.  (a)-(c) 
Real-space  particle  trajectories  at  selected  times  during  the  scan,  beginning 
(a),  middle  (b),  and  end  (c),  respectively.  Note  that  the  bottom  layer  of  the 
tip  adheres  to  the  substrate  surface  (c);  (d)  Center -of-mass  height  of  the 
tip-holder  assembly  during  the  scan,  as  a  function  of  scan  distance, 
a  «*  2.0931  A.  Note  the  decrease  in  height  associated  with  the  adherence  of 
the  bottom  tip  atoms  to  the  substrate  [see  (c)  j .  (e),  (0  Normal  force  F, 
and  tangential  force  in  the  direction  of  the  scan  F ,  during  the  scan,  (g) 
Potential  energy  of  the  tip  atoms  Ep  during  the  scan. 


son  of  the  atomic  configurations  at  the  beginning  [Fig. 
7(a)],  during  the  scan  [Fig.  7(b)  corresponding  to X—  la, 
demonstrating  a  slip]  and  towards  the  end  of  the  scan  [Fig. 
7(c)  ]  shows  that  as  a  result  of  the  interactions  between  the 
tip  and  substrate  atoms,  the  bottom  layer  of  the  tip  adheres 
to  the  substrate  and  thus  in  order  to  maintain  the  same  force 
on  the  tip  holder  throughout  the  scan  (as  required  in  the 
constant-force  scan  mode)  the  tip  assembly  must  move  clos¬ 
er  to  the  substrate.  These  simulations  demonstrate  that  in 
reactive  tip-substrate  systems,  even  under  relatively  small 
loads,  rather  drastic  structural  modifications  may  occur, 
such  as  “coating”  of  the  substrate  by  the  tip  (or  v<ce  versa). 

The  frictional  force  obtained  in  simulations  employing  a 
disordered  rigid  102-atom  tip,  prepared  by  quenching  of  a 
molten  droplet,  scanning  under  a  load  F^ex,  —  1.0  are  shown 
in  Fig.  5(e).  The  significance  of  this  result  lies  in  the  periodic 
variation  of  the  force,  reflecting  the  atomic  structure  of  the 
substrate.  This  demonstrates  that  microscopic  investiga¬ 
tions  of  structural  characteristics  and  tribological  properties 
of  crystalline  substrates  are  not  limited  to  ordered  tips.7 

Finally,  we  note  two  most  recent  studies  of  tip-substrate 
interactions  involving  other  materials.  Employing  an  em¬ 
pirically  constructed  graphite  potential  and  a  potential  for 
diamond24  (both  based  on  a  reparametrization  of  the  SW 
potentials,  and  in  the  case  of  graphite  augmented  by  inter¬ 
layer  Lennard-Jones  interactions),  it  was  shown24*  that  for 
forces  stated  in  the  experimental  literature,  a  monoatomic 
tip  would  cause  very  large  relaxations  of  the  graphite  sub¬ 
strate  extending  over  a  large  spatial  extent,®  with  an  eventual 
penetration  of  the  surface  at  5  X  10“*  N.  Furthermore,  it 
was  shown  that  by  assuming  that  the  AFM  tip  carries  a 
graphite  flake,  a  wide  variety  of  images  seen  by  the  atomic 
force  microscope  can  be  reproduced  and  correlated  with  rel¬ 
ative  rotations  of  the  flake  with  respect  to  the  scanned  sub¬ 
strate.  In  other  studies,23  employing  the  embedded  atom 
method  (EAM),  interactions  between  metal  tips  and  metal 
substrates  have  been  investigated.  These  studies  demon¬ 
strate  that  as  a  consequence  of  the  metallic  cohesion,  metal¬ 
lic  tips  lowered  onto  metallic  substrates  adhere  to  the  sub¬ 
strate  ("pressure  welding”)  resulting  in  plastic  deformation 
upon  separation  signified  by  a  pronounced  hysterises  in  the 
force  versus  distance  records  ( upon  loading  and  unloading ) , 
similar  to  the  observations  in  recent  AFM  experiments.26 


III.  DYNAMICS  OF  BOUNDED  THIN  LIQUID  FILMS 

Studies  of  nonequilibrium  systems,  flow  in  particular, 
present  immense  theoretical  and  technical  challenges.  Cou¬ 
pled  with  the  basic  conceptual  problems  presented  by  non¬ 
equilibrium  phenomena  is  the  technological  motivation  re¬ 
lated  to  the  molecular  design  of  fluids  of  desired 
rheological27  and  hydrodynamical  characteristics  for  lubri¬ 
cation  purposes.  Lubricating  fluids  are  made  usually  of  poly¬ 
meric  molecules  (sometime  in  colloidal  suspension)  which 
are  characterized  by  structural  and  shear-stress  relaxation 
times  ranging  from  milliseconds  to  minutes,  and  the  relevant 
flow  rates  are  of  the  order  of  the  reciprocal  of  these  times. 
Thus,  the  Weissenburg  number27  for  these  materials  (the 
product  of  shear  rate  and  relaxation  time)  is  of  the  order  of 
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unity.  While  simulations  of  phenomena  on  the  above  time 
scales  are  computationally  prohibitive,  it  is  possible  to  study 
the  nature  of  the  rheological  processes  occurring  in  these 
systems  via  simulations  of  simpler  fluids  sheared  at  the  same 
Weissenburg  number.2*  Simple  Lennard-Jones  systems  ex¬ 
hibit  such  “corresponding  state”  behavior. 

To  investigate  the  nature  and  properties  of  bounded  thin 
fluid  films  at  equilibrium  and  under  shear,  we  have  conduct¬ 
ed  systematic  MD  simulations  of  such  systems  under  var¬ 
ious  external  conditions  and  for  a  variety  of  interaction  pa¬ 
rameters.  We  demonstrate  our  results  for  a  system 
composed  of:  (i)  two  fee  solid  slabs  each  consisting  of  three 
dynamic  and  three  static  layers  (with  70  particles/layer) 
exposing  the  ( 1 1 1 )  surface  with  the  particles  interacting  via 
a 6-12  Lennard-Jones  (U)  potential  [seeEq.  (1),  Sec.  IV]. 
While  periodic  boundary  conditions  are  imposed  in  direc¬ 
tions  parallel  to  the  surface  plane  the  system  (dynamic  solids 
and  confined  fluid)  thickness  in  the  [  1 1 1  ]  direction  evolves 
dynamically  under  a  load  of  0.2S  e a/cr‘a  (where  ea  and  a„ 
are  the  LJ  parameters  for  the  intersolid  interactions),  and 
(ii)  ISO  fluid  particles  of  equal  size,  bounded  between  the 
dynamic  solid  interfaces,  and  interacting  between  them¬ 
selves  via  a  U  potential  with  a  well  depth  parameter  eu  =  eM 
/10,  and  with  the  solid  by  eu  =  (e„e„ ) ,/2.  The  temperature 
of  the  system  is  controlled  at  T=  0.074  (in  reduced  LJ  units, 
i.e.,  £a ),  slightly  above  the  melting  temperature  of  the  fluid, 
via  scaling  of  particles’  velocities  in  the  two  dynamic  layers 
which  are  not  in  direct  contact  with  the  liquid  in  each  of  the 
bounding  solid  slabs,  which  thus  serve  as  a  heat  bath  (i.e., 
particles  in  the  solid  layers  at  the  interface  and  the  fluid 
particles  evolve  freely  with  no  temperature  control,  transfer¬ 
ring  heat  to  the  substrate  slab  in  a  natural  manner). 

The  density  profile  in  the  direction  normal  to  the  interface 
( [  1 1 1  ],  Z  direction)  at  equilibrium  ( T  —  0.074)  isshown  in 
Fig.  8(b)  exhibiting  layering  in  the  fluid,29  and  short-time 
real-space  trajectories  of  the  solid  and  fluid  layers  at  the  in¬ 
terface  exhibiting  intralayer  order  and  interlayer  registry  are 
shown  in  Fig.  9(a).  The  slab  thickness  (in  units  ofl<rm )  vs 
time  (in  reduced  LJ  time  units,  tu)  is  shown  in  Fig.  8(a) 
with  the  end  of  the  equilibrium  simulation  corresponding  to 
t  <  75.  At  time  t  =  75  we  apply  a  sudden  shear  velocity  of  0. 1 
<7„/tu  along  the  f  1,1,0]  (x)  direction  to  the  static  portions 
of  the  solid  slabs  (which  as  a  result  start  translating  in  oppo¬ 
site  directions  since  the  velocities  are  applied  to  the  opposing 
solid  slabs  in  opposite  senses).  As  seen  from  Fig.  8(a),  upon 
application  of  the  shear  perturbation,  the  thickness  of  the 
system  initially  increases  slightly  (in  major  part  in  the  fluid 
region  which  is  set  in  motion  in  response  to  the  shear  pertur¬ 
bation)  followed  first  by  a  steady  decrease  and  subsequently 
by  an  abrupt  drop  (at  f»165).  Inspection  of  the  density 
profiles  at  an  intermediate  time  — r  =  160  [Fig.  8(c)]  and 
after  steady-state  is  achieved  [fs210.  Fig.  8(d)]  reveals 
that  the  change  in  the  fluid  film  thickness  is  due  to  a  collapse 
of  the  system  from  three  to  two  fluid  layers,  which  corre¬ 
sponds  to  an  increase  in  areal  density  of  the  fluid  layers  by 
—  20%.  The  pronounced  intralayer  structure  of  the  fluid 
during  steady-state  shear  is  shown  in  Fig.  9(b).  Accompa¬ 
nying  this  structural  transition  in  the  fluid  is  a  marked  in¬ 
crease  in  the  peak-to-peak  fluctuations  of  t  he  internal  shear 


Fig.  8.  Results  of  MD  simulations  of  a  fluid  confined  between  parallel  solid 
regions,  (a)  System  thickness  (in  units  of  7 <r„ )  vs  time  [in  LJ  tu].  For 
t  <  70  system  at  equilibrium  at  T*  0.074  f„ .  For  r  >  75  system  under  shear, 
(b)  Density  profile  in  the  z  ([111])  direction  of  the  sheared  system  at 
r  k  130  showing  three  dynamic  solid  layers  on  each  side  of  a  layered  fluid 
(middle  three  layers).  The  middle  fluid  layer  contains  —20  particles,  (c) 
Same  as  (b)  but  under  shear  at  / »  160  [see  (a)],  (d)  Same  as  (a)  but  at 
steady  state  under  shear  (/  —  210)  exhibiting  the  transition  to  a  two-layer, 
structured  fluid,  t  distance  in  (b)-(d)  in  units  of  <r„.  (e)  Internal  shear 
stress  ( a„ ,  in  units  of  em  /au  1 )  in  the  system  vs  time  exhibiting  an  increase 
in  peak-to-peak  fluctuation  upon  the  structural  transition  in  the  fluid. 


stress  {(rxt)  of  the  system  (calculated  from  the  instanta¬ 
neous  trajectories)  as  seen  in  Fig.  8(e).  The  very  regular 
fluctuations  in  thickness  (after  collapse)  seen  in  Fig.  8(a) 
coincide  in  frequency,  with  the  passing  of  ( 1 ,1,0)  rows  in  the 
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Fig.  9.  (a)  Short-time  real-space  trajectories  of  particles  in  the  interfacial 
dynamic  solid  (filled  circles)  and  fluid  (empty  circles)  layers,  exhibiting 
intralayer  order  and  interlayer  registry  at  equilibrium,  (b)  Same  as  (a)  at 
steady  state  under  shear.  Note  that  the  fluid  layer  under  these  conditions  is 
denser  than  at  equilibrium  [compare  to  (a)  ] .  The  positions  of  the  atoms  at 
the  start  of  the  short-time  interval  are  in  the  middle  of  the  circles. 


slip  plane,  indicating  persistant  order  on  both  sides  of  the  slip 
plane.  These  phenomena  have  similar  characteristics  to 
those  observed  recently  by  Israelachvili  using  an  extension 
of  the  surface  force  apparatus,3  indicating  a  structural  dy¬ 
namical  mechanism  which  (in  addition  to  film  drainage) 
can  be  operative  in  such  systems.  We  should  note  that  for  the 
above  shear  velocity  the  collapse  of  the  film  does  not  occur 
for  a  fluid  film  consisting  of  1 60  particles.  The  properties  and 
response  of  the  system  (including  film  drainage)  for  various 
film  thicknesses,  loads  and  shear  velocities  are  under  investi¬ 
gation.  Finally,  we  note  that  in  simulations  of  thick  fluid 
films  (e.g.,  three  times  thicker  than  in  the  above  example), 
under  low  shear  rates  we  have  found17  similar  interfacial 
layering  (involving  2-3  fluid  layers)  which  decays  into  the 
middle  of  the  fluid  film,  with  the  latter  region  exhibiting 
properties  (viscosity  and  structure)  similar  to  a  bulk  fluid 
under  similar  conditions.  However,  at  higher  shear  rates  the 
interfacial  liquid-layering  is  much  less  pronounced  and  in 
addition  non-Newtonian  flow  occurs. 

Finally,  recent  theoretical  studies  in  our  laboratory30, 
have  shown  evidence  for  sheaf -induced  ordering  in  atomic 
and  molecular  liquids  which  have  been  subsequently  ob¬ 
served  using  the  surface  force  apparatus.30*1 

These  results  and  studies  of  the  dependencies  of  the  flow 
characteristics  on  the  strength  of  the  fluid-solid  interaction 
and  the  composition  and  complexity  of  the  fluid  demon¬ 
strate  the  potential  utility  of  MD  simulations  in  tribological 
molecular  design  problems. 


IV.  DYNAMICAL  SIMULATIONS  OF  STRESS,  STRAIN 
AND  FINITE  DEFORMATIONS 

Common  observations  related  to  tribological  phenomena, 
and  of  thermomechanical  properties  and  response  of  materi¬ 
als  in  general,  are  usually  made  at  the  continuum  level.31 
Consequently,  (and  naturally)  the  development  of  theoreti¬ 
cal  understanding  of  these  phenomena  followed  the  "contin¬ 
uum  modelling”  approach.31  The  methodology  of  the  devel¬ 
opment  of  these  models  is  based  on  the  principles  of  mass, 
momentum  and  energy  balance,  and  the  formulation  of  con¬ 
stitutive  equations.  While  the  mathematical  fori  ulation  of 


classical  models  of  the  mechanical  response  of  matter  (such 
as  the  classical  theories  of  elasticity  and  hydrodynamics) 
achieved  a  high  degree  of  sophistication,  current  focus  is  on 
the  incorporation  of  knowledge  about  the  microscopic  be¬ 
havior  of  materials  in  continuum  models  which  attempt  to 
describe  macroscopic  observations.  This  is  done  via  the  in¬ 
troduction,  into  the  continuum  models,  of  a  set  of  state  vari¬ 
ables  which  provide  averaged  (coarsened)  representations 
of  the  relevant  microscopic  quantities.  In  addition  we  should 
note  that  most  applied  models  of  mechanical  response  are 
limited  to  the  elastic  (small  spatial  deformation)  and  linear 
(small  rates  of  application  of  the  external  perturbation)  re¬ 
gimes.  Attempts  to  incorporate  inelastic  response  and  non¬ 
linear  effects  result  in  a  great  ( often  prohibitive )  complexity. 

To  gain  knowledge  about  the  energetics  and  atomic  scale 
mechanisms  and  to  identify  and  determine  characteristic 
material  and  ambient  parameters  which  underlie  and  govern 
the  properties  and  response  of  materials  requires  experimen¬ 
tal  and  theoretical  methods  which  probe  physical  systems 
with  refined  spatial  and  temporal  resolution.  As  we  dis¬ 
cussed  above  molecular  dynamics  (MD)  simulations  afford 
such  investigations. 

Traditionally,  MD  simulations  have  been  employed  in 
studies  of  systems  of  fixed  shape  and  size  of  the  periodically 
replicated  c&lculational  cell  (Le.,  constant  volume  simula¬ 
tions).  More  recently,  methods  for  simulating  systems  in 
which  the  volume  and  shape  of  the  calculations!  cell  may 
vary  dynamically  have  been  developed,32-3*  opening  the  way 
to  investigations  of  a  large  number  of  materials  phenomena 
in  which  the  dynamical  freedom  of  the  system  to  change 
volume  and/or  structure  (or  phase)  is  essential.  In  addition, 
a  number  of  methods  have  been  developed  for  simulations  of 
fiow  and  hydrodynamics]  systems39-42  which  allow  detailed 
investigations  of  these  nonequilibrium  phenomena.  Using 
molecular  dynamics  techniques  various  studies  of  the  me¬ 
chanical  properties  of  solids  and  fluids  have  been  reported. 
Among  the  investigations  of  solid  systems  we  note  studies  of 
stressed  solids  and  crack  propagation/3-4*  dislocation  ener¬ 
getics  and  dynamics,46  structural  transformations  in  crystal 
lattices  under  uniaxial  tension  or  compression,3337,49,50  slid¬ 
ing  and  migration  of  grain  boundaries,37  simulations  of  plas¬ 
tic  deformations3*  and  shock  wave  dynamics,51,52  studies  of 
stressed  interfaces,3*’40  and  calculations  of  elastic  coeffi¬ 
cients  using  MD  simulations.53 

We  have  investigated3*-40  using  MD  simulations  the  mi¬ 
croscopic  dynamical  response,  deformation,  and  stress  relief 
mechanisms  at  crystalline  solid  interfaces  subject  to  exter¬ 
nally  applied  perturbations.  (Constant  thermodynamic  ten¬ 
sion,  constant  strain-rate  and  isoextemal-stress  simulations 
have  been  performed. )  The  objectives  of  these  studies  were: 
(i)  identification  of  the  mechanisms  for  solid  interfacial  sys¬ 
tems  of  deformation,  stress  accumulation  and  relief  and  the 
dynamical  response  to  external  perturbations,  (ii)  identifi¬ 
cation  of  the  dependence  of  the  above  phenomena  on  materi¬ 
al  characteristics,  such  as  bonding  strength,  atomic  sizes  and 
interface  crystallography  and  on  ambient  conditions  (ther¬ 
mally  adiabatic  versus  isothermal,  and  (iii)  the  develop¬ 
ment  and  critical  assessment  of  MD  simulation  methods  for 
investigations  of  the  above  phenomena. 
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Fig.  10.  A  schematic  of  the  calculation*!  cell.  is  the  number  oflayers  in 

the/I  (hard)  material  and  SL,  the  number  oflayers  in  the  5  (soft)  materi¬ 
al.  The  interface  is  between  layers  N^,  and  +  1.  The  direction  of  the 
applied  load  and  shear  stresses  are  indicated.  Three-dimensional  periodic 
boundary  conditions  are  employed  in  the  simulations. 


In  these  studies  we  focus  on  the  mechanisms  and  dynam¬ 
ics  of  the  response  of  a  system  containing  an  interface 
between  two  materials  to  external  stresses.  The  model  sys¬ 
tem  which  we  employ  in  our  simulations  consists  of  N  parti¬ 
cles  (Na  of  type  A  and  NB  of  type  B,  NA  +  NB  =  N)  inter¬ 
acting  via  pairwise  6-12  Lennard-Jones  (LJ)  potentials: 

>'(r)=4ea<,[[-^.]U-[-^-]6],(a^)=^,  (1) 

where  A  and  B  represent  two  types  of  materials.  The  solid  is 
set  up  initially  in  a  fee  crystalline  structure,  with  NL  (111) 
layers  (N/NL  particles  per  layer)  with  the  Z  axis  along  the 
[111]  direction,34  and  three-dimensional  periodic  boundary 
conditions  are  used  (see  Fig.  10) .  The  well-depth  parameter 
)  of  the  interaction  potential  between  particles  in  layers 
1  through  Nm  is  taken  to  be  twice  that  for  particles  in  layers 
+  1  through  N  (i.e.,  eAA  =  2e„)  corresponding  to  a 
soft,  solid,  lubricating  material  (the  B  system)  pressed 
between  hard-material  slabs  (the ,4  system).  In  order  to  iso¬ 
late  the  dependence  of  the  system  response  on  the  interaction 
strength  (potential  depth)  parameter  e  from  that  due  to  dif¬ 
ferences  in  the  atomic  sizes  ( and  thus  interatomic  distances ) 
associated  with  the  parameter  a  in  the  LJ  potential,  we  have 
performed  first  simulations  in  which  eAA  =  2eBg  and  aAA 
—  crBB,  and  then  simulations  in  which  aBB  —  l.Scr^.  The 
interspecies  LJ  potential  parameters  were  chosen  as  aAB 
=  {<*aa  +  crBB)/2zn&€AB  =  (eAAeBB ) ,/J.  In  the  following 
we  use  reduced  units33  where  energy  and  temperature  are 
expressed  in  units  of  e M ,  length  in  units  of  aAA ,  stress  in 
units  of  (cAA/<?AA )  and  the  time  unit  (tu)  is  (mA/eAA  )l/2 
aAA .  In  the  integration  of  the  equations  of  motion  we  used  a 
time  step  Af  =  0.0075  tu  which  results  in  energy  conserva¬ 
tion  (to  six  significant  figures)  for  extended  runs. 

Following  equilibration  of  the  initial  system  at  a  reduced 
temperature  of  0.1 1  (note  that  a  pure  bulk  LJ  crystal 
melts  at  T=*  0.7,  and  thus,  since  eBB  —  eAA  /2,  the  bulk  melt¬ 
ing  point  of  the  soft  material  B  is  half  that  of  the  A  material ) 
under  a  load  of  0.5  in  the  [  1 1 1  ]  (Z)  direction  we  apply  the 
external  perturbation.  In  simulations  of  applied  thermody¬ 
namic  tension  [isoenthalpic,  isotension  and  constant  num¬ 
ber  of  particle,  N,  simulations,  i.e.,  the  {%*,r,N)  ensemble] 
we  apply  to  the  system  a  thermodynamic  tension  tensor  r 
(which  is  related  to  the  external  stress  tensor)  in  the  [  1,1,0] 
direction  (see  Fig.  10)  at  a  rate  of  r„  =  0.00125  ( e/ai)/At 


and  follow  the  evolution  of  the  system  until  it  fails  [to  keep 
the  system  from  a  rigid  rotation  a  symmetric  thermodynam¬ 
ic  tension  tensor  r  is  applied,  (i.e.,  =  r„ )  ] . 

To  investigate  the  dependence  on  the  thermal  ambient 
conditions  we  distinguish  between  simulations  where  the 
system  is  thermally  isolated  during  the  application  of  the 
shear  and  simulations  where  isothermal  conditions  are 
maintained.  As  the  system  envolves  under  the  applied  shear 
it  develops  internal  stresses  which  are  calculated  using  the 
positions  and  forces  on  the  particles.38  These  simulations 
allow  us  to  determine  the  nature  of  structural  transforma¬ 
tions  and  the  critical  external  stresses  (or  critical  strains  in 
constant  strain-rate  simulations)  which  characterize  the 
transition  from  elastic  to  inelastic  response  regimes,  and  the 
critical  yield  stress  (where  the  internal  stresses  vanish). 

The  mechanical  response  and  energetic  characteristics  of 
the  system  can  be  investigated  via  layer  decomposition  of  the 
system  properties.  Results  at  a  constant  external  =  0.8 1 
and  adiabatic  thermal  conditions  for  a  system  composed  of 
nine  layers  each  of  equal  size  particles  of  hard  and  soft  mate¬ 
rials  (with  70  particles/layer  and  the  interface  located 
between  layers  9  and  10),  which  exhibited  yield  at  a  critical 
value  of  r,  „  =  0.86,  are  shown  in  Figs.  1 1  and  12.  Inspec¬ 
tion  of  the  real  space  trajectories  revealed  that  at  this  value  of 
stress  (which  is  below  the  critical  yield  value)  the  system 
undergoes  a  structural  transformation  involving  interlayer 
motion  leading  to  stacking  fault  generation  (between  layers 
12  and  13)  inside  the  soft  material  and  slip. 

In  Fig.  11(a)  and  1 1  (b)  we  depict  the  per  layer  xz  compo¬ 
nent  of  the  internal  stress  for  the  interface  layers  (layer  5, 8, 
and  9  of  the  hard  material  and  layers  10-14  of  the  soft  mate- 


Fio.  1 1.  Per-layer  internal  stresses  (a'„)  vs  time  (in  tu)  for  a  simulation  at 
ra  m  0.8 1  under  adiabatic  condition*.  Layers  S,  S,  and  9  in  the  hard  materi¬ 
als  and  10  and  11  in  the  soft  ones  are  shown  in  (a).  The  stresses  in  layers  12- 
14  in  the  toft  material  are  shown  in  (b).  The  interface  is  between  layer*  9 
and  10.  Note  the  variation  in  the  internal  stresses  at  the  time  of  the  struc¬ 
tural  transformation.  The  internal  stresses  in  the  interior  of  the  toft  material 
(10-14)  decrease  with  layer  10  exhibiting  pinning  by  the  hard  material.  The 
stress  relief  in  the  soft  material  it  accompanied  by  stress  accumulation  in  the 
hard  material  (layers  5-9). 
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Fig.  12.  Per  layer  potential  en¬ 
ergies,  £'in  (a),  and  tempera¬ 
tures  V  in  (b),  for  the  system 
simulated  at  =.  0.81  under 
adiabatic  conditions. 


rial).  From  these  figures  we  observe  that  the  generation  of 
the  structural  transition  involves  a  gradual  decrease  in  the 
internal  <r„  component  in  layers  1 1-14  of  the  soft  material 
while  the  variation  in  the  stress  in  layer  10  [  the  atomic  plane 
of  the  soft  material  (B)  adjacent  to  the  hard  material  (A), 
see  Fig.  10]  is  smaller.  As  seen  from  Fig.  11(a)  the  stress 
relief  process  and  the  associated  structural  transformation  in 
the  soft  material  are  accompanied  by  a  stress  accumulation  in 
the  interfacial  region  of  the  hard  material  ( layer  9 ) .  In  addi¬ 
tion  we  observe  periodic  oscillations  in  the  internal  stress 
(most  pronounced  for  layers  11-14)  past  the  structural 
transformation  period,  as  the  system  relaxes  in  the  new  state 
after  the  structural  transformation  events. 

The  energetics  of  the  system  is  explored  via  the  time  re¬ 
cords  of  the  perlayer  potential  energies  (£J)  and  tempera¬ 
ture  (7*)  shown  in  Figs.  12(a)  and  12(b),  respectively. 
From  the  potential  energy  curves  we  find  that  the  potential 
energy  of  particles  in  layer  1 1  (second  layer  from  the  inter¬ 
face,  see  Fig.  10)  in  the  soft  material  is  initially  lower  than 
that  of  layers  12-14,  which  are  further  removed  from  the 
interface,  since  particles  in  that  layer  are  within  the  range  of 
the  stronger  interaction  with  the  hard  substrate 
[  fAM  *  *aa  *mb  ) 1/2  ]  •  The  potential  energy  of  particles  in  the 
interfacial  soft  material  layer  (layer  10)  adjacent  to  the  hard 
material  is  lower  by  about  —  0.4  e  than  the  value  for  layer 
1 1.  This  extra  stabilization  “pins”  this  interfacial  layer  to  the 
hard  substrate.  The  potential  energy  of  the  topmost  layer  of 
the  hard  material  ( layer  9 )  is  lower  further  by  about  —  2.6e. 
The  potential  energies  of  both  layers  9  and  10  increase  due  to 


the  structural  transformation.  Since  the  system  in  this  set  of 
calculations  is  thermally  isolated,  the  structural  change  is 
accompanied  by  a  temperature  increase  as  seen  from  Fig. 
12(b)  (where  the  curves  for  successive  layers,  starting  from 
10  and  up,  are  displaced  vertically  by  0.04e).  Note  however 
that  the  final  temperature  after  the  transformation  is  below 
the  melting  temperature  for  both  materials. 

To  investigate  the  dependence  of  the  system  properties  on 
the  mismatch  in  atomic  sizes  between  the  two  interfacing 
materials,  we  have  performed  simulations  in  which  the  well- 
depth  parameters  were  chosen  as  before  but  the  o  param¬ 
eters  were  chosen  such  that  the  atoms  in  the  soft  ( B )  materi¬ 
al  are  characterized  as  having  a  larger  size,  i.e.,  aBB  =  1.5 
oAA  and  oAB  —  1.25  oAA.  From  these  simulations  we  have 
determined  that  the  critical  yield  value  of  the  thermodynam¬ 
ic  tension  in  this  system  is  rcja  =  0.62  ±  0.01  (and  the  cor¬ 
responding  external  stress  crKja  =  0.51  ±  0.01 ),  which  are 
considerably  lower  than  the  corresponding  values  found  for 
the  equal-atomic  size  systems.  Moreover,  unlike  the  pre¬ 
vious  cases,  where  yield  was  preceded  by  an  external  stress 
regime  in  which  the  system  responded  inelastically  via  a  se¬ 
quence  of  structural  transformations,  for  the  present  system 
(in  which  the  atomic  sizes  across  the  interface  differ)  no 
such  behavior  is  observed.  Furthermore,  the  whole  soft  sys¬ 
tem  responds  in  unison,  starting  at  the  interfacial  layer  and 
up  into  the  material.  These  observations  can  be  understood 
when  considering  that  as  a  consequence  of  the  larger  atomic 
size  the  atoms  of  the  soft  material  at  the  interface  average 
over  the  corrugation  of  the  potential  due  to  the  substrate, 
resulting  in  an  effective  potential  surface  which  exhibits . 
smaller  variations  for  lateral  displacements  parallel  to  the 
interface  plane,  and  consequently  a  reduced  resistance  to 
shear. 

The  results  of  our  extensive  investigations3*  may  be  sum¬ 
marized  as  follows:  (i)  For  interfacial  systems  which  are 
characterized  by  differing  interatomic  interaction  strengths 
(i.e.,  the  interface  is  between  a  hard  and  soft  material),  the 
system  responds  to  an  applied  nonisotropic  perturbation 
(applied  shear  or  strain)  first  elastically  and  then  via  stress 
relief  mechanisms  which  involve  structural  transformations 
(stacking  fault  formation  and  interlayer  slip).  For  larger 
values  of  the  external  forces,  eventual  yield  occurs,  (ii)  Criti¬ 
cal  values  of  the  external  perturbation  required  in  order  to 
bring  about  inelastic  response  (structural  transformations 
and  yield)  have  been  determined.  Our  simulations  demon¬ 
strate  that  the  critical  values  are  smaller  for  a  system  under 
thermally  adiabatic  conditions  than  for  an  isothermal  envi¬ 
ronment  The  origin  of  the  difference  lies  in  the  dissipation  of 
the  generated  thermal  energy  under  isothermal  conditions, 
which  necessitates  larger  values  of  the  external  perturba¬ 
tions  in  order  to  overcome  potential  barriers  for  structural 
modifications  and  eventual  yield,  (iii)  The  cohesive  intera¬ 
tomic  interactions  at  the  interface  between  a  hard  substrate 
and  a  soft  material  result  in  “pinning”  of  the  soft  material  at 
the  interface  ( 1-3  atomic  layers) .  As  a  result  the  response  of 
the  system  to  the  external  perturbation  ( stress  relief  via  plas¬ 
tic  as  well  as  structural  transformations)  occurs  in  a  “shear 
band”  consisting  of  a  few  atomic  layers  inside  the  soft  mate¬ 
rial,  which  for  our  model  system  are  located  at  about  1-3 
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layers  away  from  the  original  (unstressed)  interface.  The 
stress  relief  in  the  soft  material  is  accompanied  by  stress  ac¬ 
cumulation  in  the  hard  substrate,  (iv)  MD  simulations  for 
interfacing  hard  and  soft  materials,  which  in  addition  are 
characterized  by  differing  atomic  sizes,  reveal  the  important 
role  played  by  atomic  size  mismatch  in  determining  the 
atomic-scale  mechanism  of  response.  For  such  system  it  was 
found  that  no  adhesive  “pinning”  occurs  at  the  interface, 
and  that  the  soft  (and  larger  atomic  size)  material  responds 
as  a  whole  with  no  distinct  structural  transformations  pre¬ 
ceding  the  yield  point.  The  critical  yield  stress  value  for  this 
system  is  significantly  lower  than  that  found  for  the  corre¬ 
sponding  equal-atomic-size  system,  (v)  Comparison  of  our 
results  in  this  study  for  the  [  1 1 1  ]  interface  with  our  previous 
investigations  of  the  [  100]  interface40  demonstrates  the  de¬ 
pendence  of  the  critical  values  of  the  shear  stresses  on  the 
crystallographic  orientation  of  the  interface,  as  well  as  of 
certain  details  of  the  stress  relief  mechanisms. 
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Dynamical  interactions  between  a  scanning  tip  and  a  silicon  substrate  are  investigated  using 
molecular  dynamics  simulations  of  both  the  constant-height  and  constant-force  scan  modes. 
Localized  temporary  and  permanent  modifications  of  the  substrate  occur,  depending  on  tip-sub¬ 
strate  separation  and  scan  geometry.  Implications  for  resolving  structural  and  force  characteristics 
in  scanning  tip  spectroscopies,  employing  atomically  sharp  as  well  as  large  ordered  of  disordered 
tips  are  discussed. 


The  developments  of  scanning  tunneling  microscopy  [lj  (STM),  the  related 
atomic  force  microscopy  (AFM)  [2],  and  the  surface  force  apparatus  (SFA) 
[3,4]  revolutionized  our  perspectives  and  abilities  to  probe  the  morphological 
and  electronic  structure  and  the  nature  of  interatomic  forces  in  materials,  as 
well  as  opened  new  avenues  [5]  for  microscopic  investigations  and  manipula¬ 
tions  of  technological  systems  and  phenomena  such  as  tribology  [4.6],  lithogra¬ 
phy  and  in  biochemical  applications.  In  both  STM  and  AFM  a  tip  is  brought 
close  to  the  surface  and  either  the  tunneling  current  (STM)  or  deflection  of  the 
cantilever  holding  the  tip  (AFM)  are  monitored  as  the  surface  is  scanned.  Of 
particular  interest  are  questions  related  to  the  nature  of  the  tip-substrate 
interactions  and  the  dynamical  response  of  the  substrate  (and  tip)  which  may 
result  in  temporary  or  permanent  modifications  of  the  local  properties  [7,8] 
and  be  reflected  in  the  recorded  data. 

In  this  paper  we  investigate  the  dynamics  of  the  tip- substrate  interactions 
via  molecular  dynamics  simulations  employing  a  realistic  interatomic  interac¬ 
tion  potential  for  silicon  (the  Stillinger-Weber  (SW)  potential  [9])  which  has 
been  used  recently  in  a  number  of  investigations  [10-14]  of  bulk,  surface  and 
interface  properties  of  this  material.  Since  both  the  substrate  and  tip  consist  of 
the.  same  material,  these  simulations  correspond  to  the  case  of  a  reactive 
tip-substrate  system  [11-14],  Our  simulations  [12-14],  in  both  the  constant- 

*  Work  supported  by  US  DOE  under  grant  No.  FG05-86ER45234.  and  by  the  DARPA-Hughes 
Tribology  Program. 
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tip-height  and  constant-force  scan  modes,  reveal  that  the  local  structure  of  the 
surface  can  be  stressed  and  modified  as  a  consequence  of  the  tip-substrate 
dynamical  interaction,  even  at  tip-substrate  separations  which  correspond  to 
weak  interaction.  For  large  separations  these  perturbations  anneal  upon 
advancement  of  the  tip  while  permanent  damage  can  occur  for  smaller 
separations.  For  the  material  that  we  simulated  (Si),  we  do  not  find  long-range 
elastic  deformations,  which  may  occur  in  other  circumstances  [7]  depending 
upon  the  elastic  properties  of  the  material  and  the  nature  of  interactions.  The 
characteristics  of  the  data  depend  upon  the  geometry  of  the  scan,  the  degree  of 
perfection  of  the  substrate  and  the  temperature.  We  identify  various  dynami¬ 
cal  events  including  stick-slip  phenomena,  which  could  be  experimentally 
resolved,  using  current  estimates  [1,2,5],  and  which  would  influence  the 
analysis  of  data,  as  well  as  pointing  to  ways  in  which  the  tip  could  be  used  for 
atomic  scale  manipulations  of  the  material. 

In  our  simulations  the  system  consists  of  4  layers  of  dynamic  Si  particles 
with  49  (or  100)  atoms  per  layer,  exposing  the  (111)  surface,  and  interacting 
with  2  layers  of  a  static  Si  substrate  of  the  same  structure  (calculations  with 
systems  of  6  dynamic  layers  and  a  larger  number  of  particles/layer  yield  very 
similar  results).  The  2D  calculational  cell,  defined  by  the  (110)  and  (101) 
vectors,  is  periodically  repeated  parallel  to  the  (111)  plane.  Two  types  of  tips 
were  employed  in  our  simulations.  First,  we  have  used  a  sharp  tip  [1,2,5] 
simulated  by  4  Si  atoms  in  an  initial  tetrahedral  configuration,  “mounted”  on 
and  interacting  with  2  layers  of  silicon  atoms  which  serve  as  a  holder.  In 
addition  we  have  employed  larger  tips  consisting  of  102  atoms,  which  were 
either  ordered  (exposing  a  16  atoms  (111)  planar  facet)  or  disordered.  The 
equations  of  motion,  governed  by  the  SW  potential  [9],  which  contains  2-  and 
3-body  interactions  {V2  and  V}  respectively),  are  integrated  using  a  5th  order 
predictor-corrector  algorithm  with  a  time  step  At  — 1.15  x  10-3  ps  (or  0.015 
tu  where  tu  =  7.6634  x  10  “ 14  s).  Throughout  we  use  e  =•  50  kcal/mol  as  the 
unit  of  energy,  a  —  2.0951  A  as  the  unit  of  length,  and  t/o  —  1.65728  X  10" 9 
N  as  the  unit  of  force.  The  reduced  units  of  length  X*,  Y *  and  Z*,  along  the 
(110),  (101)  and  (111)  directions  are  12.82o,  12.82a,  and  5.8333a,  respectively. 
The  kinetic  temperature  is  controlled  via  scaling  of  particle  velocities  in  the 
bottom  layer  of  the  dynamic  substrate. 

Both  constant-tip-height  and  constant-force  scan  modes  were  simulated.  In 
the  first  mode  following  equilibration  at  room  temperature,  with  the  tip 
outside  the  range  of  interaction,  the  tip  is  lowered  slowly  (5.4  x  10-4  a/Ar)  to 
a  prescribed  height.  Studies  at  3  initial  tip  heights,  A,  (/  -  1,  3),  corresponding 
to  distances  of  2.91,  2.345  and  1.227  A  between  the  lowest  tip  atom  and  the 
uppermost  layer  of  the  substrate,  were  performed.  h2,  h2  and  h2  correspond 
to  the  attractive,  equilibrium,  and  repulsive  regions  of  the  interparticle  poten¬ 
tial,  respectively.  To  faithfully  simulate  the  laboratory  process,  which  is  much 
slower  than  can  be  achieved  in  computer  simulations,  and  record  data  for 


U.  Landman  et  at.  /  Dynamics  of  tip-substrate  interactions  in  AFM 


LI  79 


structurally  relaxed  tip-substrate  configurations  we  have  adopted  the  follow¬ 
ing  scan  and  data  accumulation  procedure:  (a)  The  system  is  first  equilibrated 
for  200  A /  with  the  tip  at  the  desired  height,  (b)  Lateral  scans  (along  X* 
(110))  consist  of:  (i)  equilibration  for  150  A /  followed  by  150  Ar  of  data 
collection  and  averaging,  (ii)  motion  of  the  tip  assembly  for  50  At  to  a  new 
scan  point  (8  scan  points  per  2D  unit  cell  length,  covering  in  each  lateral  scan 
3  unit  cells),  after  which  step  (i)  is  repeated.  In  exploratory  studies  we  found 
that  the  results  of  the  simulations  are  not  modified  in  any  substantial  way  by 
increasing  the  equilibration  periods.  Results  for  2  lateral  scans  at  each  height 
are  described,  denoted  by  />,(1)  for  the  scan  on  top  of  an  atomic  row,  and  by 
>i,(2)  for  the  scan  mid-distance  between  rows. 

In  the  constant-force  simulations  in  addition  to  the  particle  equations  of 
motions  the  center  of  mass  of  the  tip-holder  assembly,  Z,  is  required  to  obey 
MZ -  (F(t)  -  Fe%t)‘Z- yZ  where  F  is  the  total  force  exerted  by  the  tip 
atoms  on  the  static  holder  at  time  t,  which  corresponds  to  the  force  acting  on 
the  tip  atoms  due  to  their  interaction  with  the  substrate,  Fat  is  the  desired 
(prescribed)  force  for  a  given  scan,  y  is  a  damping  factor  and  M  is  the  mass  of 
the  holder.  In  these  simulations  the  system  is  brought  to  equilibrium  for  a 
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Fig.  1.  (a-d)  Particle  trajectories  for  the  low-tip,  /i3(l),  scan:  (a.  d)  complete  scan;  (b.  c)  at 
beginning  and  towards  the  end  of  the  scan,  respectively.  In  (a-c)  view  along  the  (101)  direction,  in 
(d)  view  along  (111),  (e-h)  same  as  (a-d)  for  between-the-rows.  h}(2).  scan;  note  the  damage 
caused  in  the  substrate,  (i.j)  trajectories  at  the  beginning  and  towards  the  end  of  the  high-tip. 
A,(l)  scan.  Arrows  in  (c)  and  (j)  point  to  induced  defects.  The  tip  atoms  are  located  on  the 
horizontal  lines  above  the  substrate,  scanning  from  left  to  right. 
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prescribed  value  of  Fcxl,  and  the  scan  proceeds  as  described  above  while  the 
height  of  the  tip-holder  assembly  adjusts  dynamically  according  to  the  above 
feed-back  mechanism. 

In  fig.  1  we  show  particle  trajectories,  viewed  from  the  side  along  the  (101) 
direction  and  from  the  top  (along  the  (111)  direction,  d  and  h)  for  constant 
tip-height  simulations  and  a  static  tip  (i.e.,  hard  but  reactive  tip),  at  two  tip 
heights,  hx  and  /i3.  Trajectories  for  a  complete  scan  at  the  low  tip  height 
Aj(l),  are  shown  in  a  and  d  (the  scan  direction  is  from  left  to  right).  The 
horizontal  lines  are  the  static  tip-particle  trajectories  (only  3  tip  particles  out 
of  the  4  are  visible  from  this  viewing  direction).  In  figs,  lb  and  lc  particle 
trajectories  are  shown  at  the  beginning  and  towards  the  end  of  the  scan.  As 
seen  in  fig.  lb  the  interaction  with  the  tip  triggers  local  displacements  of 
substrate  particles.  In  particular,  the  atom  in  the  top  layer  right  below  the  tip 
drops  to  an  interstitial  position  in  the  second  layer.  This  is  seen  clearly  in  fig. 
lc,  where  the  interstitial  atom  is  marked  by  an  arrow.  Note  that  once  the  tip 
advances,  the  localized  defect  induced  by  the  tip  anneals  (compare  the  left 
part  of  the  scan  in  figs,  lb  and  lc).  In  comparison,  locating  the  tip  at  A,(l),  in 
the  region  of  attractive  interaction,  results  in  a  localized  outward  displacement 
of  atoms  in  the  substrate  top  layer  (see  figs,  li  and  lj,  where  the  displaced 
particle  is  marked  by  an  arrow).  The  dependence  on  the  scan  geometry  which 
affects  the  relative  tip  to  substrate  atoms  bond  distances,  orientations  and 
coordination  and  subsequently  the  resultant  forces,  is  demonstrated  by  com¬ 
paring  the  trajectories  shown  in  figs,  la-ld  and  le-lh.  As  seen,  the  between- 
row  scan,  /i  3(2),  produces  a  more  significant  damage  to  the  surface  due  to  the 
bonding  geometry  generating  top  layer  vacancies  (2  atoms  attached  to  the  tip, 
figs,  le  and  lg). 

Records  of  the  forces  on  the  tip  atoms  for  the  three  scan  heights  versus  tip 
location  and  versus  time  are  shown  in  figs.  2a-2f  and  2g  and  2h,  respectively. 
We  start  with  the  time-histories  of  the  h^l)  and  h3(  1)  scans  (figs.  2g  and  2h). 
For  the  high-tip  scan,  h^l)  (fig.  2g)  the  tip  experiences  throughout  an 
attractive  force  ( Ft  <  0)  towards  the  substrate,  increasing  as  the  tip  advances 
laterally  in  the  regions  between  substrate  atoms.  A  different  situation  is 
encountered  in  the  low-tip  scan,  A3(  1)  (fig.  2h).  Here,  as  the  tip  is  lowered  it 
initially  attracts  to  the  surface,  followed  by  a  decrease  in  the  force  as  the  tip 
enters  the  repulsive  region,  and  subsequently  experiencing  a  repulsive  interac¬ 
tion  (Fz>  0  at  ~  30  tu).  This  is  then  followed  by  a  sharp  attraction  ( -  40  tu), 
culminating  in  a  resultant  attractive  force  (Ft*  -1  e/cr),  although  the  tip  at 
this  point  is  in  intimate  contact  with  the  surface,  which  would  have  been 
expected  to  yield  a  strong  repulsion.  The  dynamical  mechanism  underlying  the 
observed  attraction  is  the  generation  of  surface  interstitial  defects  (see  also 
figs,  la-ld).  Thus,  the  interaction  with  the  tip  can  induce  local  rearrangement 
in  the  substrate  and  consequently  can  alter  (even  reverse  the  sign)  of  the 
resultant  recorded  forces.  Further  lateral  scanning  over  substrate  atoms  is 
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Fig.  2.  (a-c)  Forces  on  the  tip  atoms  for  the  high,  intermediate  and  low  tip  scans  along  the  X 
(llO)  direction  for  scan  1.  (d)  Forces  on  the  tip  for  the  between-the-rows  scan  (2);  the  F,  force  on 
the  tip  in  a  A2(l)  scan  for  a  defective  surface  is  shown  in  (e)  and  in  an  A,(l)  scan  for  a  dynamic 
<>P  (0:  (g.h)  Ft  force  on  the  tip  in  high,  A,(l),  and  low,  A}(  1),  scans  versus  time,  respectively. 
These  time  records  include  the  stage  of  lowering  the  tip.  All  others  start  after  equilibration  at  the 
indicated  height  Distance  in  units  of  a,  time  in  units  of  tu  and  force  in  units  of  e/a  - 1.6572 

xlO-9  N. 

characterized  by  periodic  oscillations  in  the  force  which  reflect  the  substrate 
periodicity. 

The  dependence  of  the  character  of  the  recorded  forces  versus  position  on 
tip  height  is  illustrated  in  figs.  2a-2c  for  scan  1.  These  records  were  obtained 
after  the  initial  equilibration,  with  the  tip  lowered  to  the  desired  height. 
Significant  variations  in  the  magnitudes  and  character  of  the  forces  are 
observed  depending  upon  the  tip  separation  from  the  surface.  For  example, 
the  Ft  component  at  h,(l)  is  overall  attractive  with  sharp  negative  spikes 
which  for  the  intermediate  tip  height,  h2(  1),  turn  into  broad  attractive 
plateaus  when  the  tip  is  scanning  between  surface  atoms.  For  the  low-tip 
configuration,  fig.  2c,  passage  between  surface  atoms  is  exhibited  by  broad 
repulsive  peaks.  Spatial  variations  of  the  potential  energy  surface  are  exhibited 
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by  comparison  of  the  two  scans  at  the  hx  height  (figs.  2a  and  2d).  Note  in 
particular  the  changes  in  F.,  and  the  change  in  magnitudes  of  Fx  and  Fy. 
These  characteristics  portraying  the  surface  structure  and  the  spatial  variation 
of  the  potential  are  within  the  current  estimates  of  resolution  in  AFM  [2], 

A  comparison  between  the  F2  records  for  an  ideal  surface  (fig.  2b)  and  for  a 
defective  surface  (fig.  2e)  reveals  the  scale  and  character  of  variations  caused 
by  microscopic  modifications  in  the  surface  morphology.  Here,  the  surface  was 
initially  equilibrated  with  a  top  layer  vacancy  at  the  location  of  the  fourth 
atom  from  the  right  in  the  row  directly  under  the  lower-most  tip  atom  (see  fig. 
Id).  Examination  reveals  that  in  course  of  the  scan  the  tip  caused  major 
rearrangement  of  the  defective  region  via  a  complicated  sequence  of  atomic 
displacements  involving  first  and  second  layer  atoms. 

Simulations  were  also  performed,  at  constant-height,  for  tip  atoms  which 
were  allowed  to  evolve  dynamically.  The  recorded  Fz  force  for  a  high-tip  scan, 
/tj(l),  shown  in  fig.  2f,  exhibits  a  decrease  in  magnitude  originating  from  a 
relaxation  of  the  tip,  occurring  mainly  at  the  beginning  of  the  scan,  and  shifts 
and  pronounced  asymmetry  of  characteristic  features  (compare  to  F:  in  fig. 
2a),  connected  with  a  stick-slip  behavior  exhibited  by  the  lower-most  tip  atom 
as  the  scan  progresses  (see  below). 

We  turn  next  to  constant-force  simulations  employing  an  initially  ordered  tip 
consisting  of  102  atoms  (exposing  a  16  atom  (111)  facet)  scanning  a  6  layer 
substrate  with  100  atoms/layer.  In  these  simulations  both  the  substrate  and 
tip  atoms  respond  dynamically  and  the  scan  rate  is  half  of  that  described 
above.  In  fig.  3  we  show  results  for  a  scan,  for  a  constant  force  value 
Fx  ex,  -  -13.0  (i.e.,  2.15  x  10-8  N).  Side  views  of  the  system  trajectories  at  the 
beginning  and  end  stages  of  the  scan  are  shown  in  figs.  3a,  3b  and  3c, 
respectively.  As  seen  the  tip-substrate  interaction  induces  local  modifications 
of  the  substrate  and  tip  structure,  which  are  transient  and  similar  to  those 
observed  in  the  constant-height  simulations.  The  recorded  force  on  the  tip- 
holder  in  the  X*  direction  is  shown  in  fig.  3d.  While  the  normal  force,  F., 
fluctuates  around  the  prescribed  value  as  required,  the  force  along  the  scan 
direction  (Fx,  in  fig.  3d)  exhibits  a  periodic  modulation  portraying  the  peri¬ 
odicity  of  the  substrate.  Most  significant  is  the  stick-slip  behavior  signified  by 
the  asymmetry  in  Fx  (observed  also  in  the  real-space  atomic  trajectories  in 
figs.  3a  and  3b).  Here,  the  tip  atoms  closest  to  the  substrate  attempt  to  remain 
in  a  favorable  bonding  environment  as  the  tip-holder  assembly  proceeds  to 
scan.  When  the  forces  on  these  atoms  due  to  the  other  tip  atoms  exceed  the 
forces  from  the  substrate,  they  move  rapidly  to  minimize  the  Fx  force,  by 
breaking  their  current  bonds  to  the  surface  and  forming  new  bonds  in  a  region 
translated  by  one  unit  cell  along  the  scan  direction.  We  note  that  our 
constant-force  simulation  method  corresponds  to  the  experiments  in  ref.  [6]  in 
the  limit  of  a  stiff  wire  (lever)  and  thus  the  stick-slip  phenomena  which  we 
observed  are  a  direct  consequence  of  the  interplay  between  the  surface  forces 
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Fi*.  3.  (a  -c)  Particle  trajectories  in  a  constant-force  simulation  viewed  along  the  (101)  direction 
just  before  (a)  and  after  (b)  a  stick-slip  event  and  towards  the  end  of  the  scan  (c),  for  a  large, 
initially  ordered,  dynamic  tip.  (d)  The  recorded  Fx,  exhibiting  stick-slip  behavior,  (e)  The  Fx 
force  in  a  constant-force  scan  (  F.  at  - 1.0)  employing  a  glassy  static  tip,  exhibiting  the  periodicity 
of  the  substrate.  Shown  in  the  inset  are  the  real-space  trajectories  towards  the  end  of  the  scan, 
demonstrating  the  tip-induced  substrate  local  modifications. 


and  interatomic  interactions  in  the  tip.  The  Fx  force  which  we  record  corre¬ 
sponds  to  the  frictional  force.  From  the  extrema  in  Fx  (fig.  3d)  and  the  load 
(Fjext)  “sed  we  obtain  a  coefficient  of  friction  p  -  |  Fx  |/|  Fznt  |  =  0.77,  in  the 
range  of  typical  values  obtained  from  tribological  measurements  in  vacuum. 

Finally,  we  show  in  fig.  3e  the  frictional  force  obtained  in  simulations 
employing  a  disordered  static  102-atom  tip,  prepared  by  quenching  of  a 
molten  droplet,  scanning  under  a  load  Fs  wl=«1.0.  The  significance  of  this 
result  lies  in  the  periodic  variation  of  the  force  reflecting  the  atomic  structure 
of  the  substrate.  This  demonstrates  that  microscopic  investigations  of  struct¬ 
ural  characteristics  and  tribological  properties  of  crystalline  substrates  are  not 
limited  to  ordered  tips  [6], 
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Abstract 

Molecular  dynamics  simulations  and  atomic  force  microscopy  are  employed 
to  investigate  the  atomistic  mechanisms  of  adhesion,  contact  formation,  nano¬ 
indentation,  separation  and  fracture  which  occur  when  a  nickel  tip  interacts 
with  a  gold  surface.  The  theoretically  predicted  and  experimentally  measured 
hysteresis  in  the  force  versus  tip-to-sample  distance  relationship,  found  upon 
approach  and  subsequent  separation  of  the  tip  from  the  sample,  is  related  to 
inelastic  deformation  of  the  sample  surface  characterized  by  adhesion  of  gold 
atoms  to  the  nickel  tip  and  formation  of  a  connective  neck  of  atoms.  At  small 
tip-sample  distances,  mechanical  instability  causes  the  tip  and  surface  to 
Jump-to-contact  which  leads  to  adhesion  induced  wetting  of  the  nickel  tip  by 
gold  atoms.  Subsequent  Indentation  of  the  substrate  results  in  the  onset  of 
plastic  deformation  of  the  gold  surface.  The  atomic-scale  mechanisms  under¬ 
lying  the  formation  and  elongation  of  a  connective  neck,  which  forms  upon 
separation,  consist  of  structural  transformations  involving  elastic  and 
yielding  stages. 
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INTRODUCTION 

Understanding  the  atomistic  mechanisms,  energetics  and  dynamics 
underlying  the  interactions  and  physical  processes  which  occur  when  two 
materials  are  brought  together  (or  separated)  is  fundamentally  important  to 
basic  and  applied  problems  such  as  adhesion  [1-7],  contact  formation  [3-16], 
surface  deformations  [17-22,7,16],  materials  elastic  and  plastic  response 
characteristics  [17-22],  materials  hardness  [23-25],  microindentation 
[6,10,24-27],  friction  and  wear  [ 16, 18b, 28-30]  and  fracture  [31,32].  These 
considerations  have  motivated  for  over  a  century  [17,18,1,3]  extensive 
theoietical  and  experimental  research  endeavors  of  the  above  phenomena  and 
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their  technological  consequences.  Most  theoretical  approaches  to  these 
problems,  with  a  few  exceptions  [7,14-I6]yhave  been  anchored  in  continuum 
elasticity  and  contact  mechanics  [17-22].  Similarly,  until  quite  recently 
[33-36]  experimental  observations  and  measurements  of  surface  forces  and  the 
consequent  materials  response  to  such  Interactions  have  been  macroscopic  in 
nature . 

The  everlasting  quest  to  understand  and  observe  natural  phenomena  on 
refined  microscopic  scales  has  led  to  the  development  of  conceptual  and 
technological  devices  allowing  the  interrogation  of  materials  with  Increasing 
resolution.  On  the  experimental  front  the  developments  of  the  surface  force 
apparatus  (SFA)  [34],  of  scanning  tunneling  microscopy  (STM)  [35]  and  of  the 
related  atomic  force  microscopy  (AFM)  [33],  broaden  our  perspectives  and 
abilities  to  probe  the  morphology,  electronic  structure  and  nature  of 
interatomic  forces  in  materials,  as  well  as  enhance  our  ability  to  manipulate 
materials  on  the  atomic  scale  [36]. 

On  the  theoretical  front,  recent  advances  in  the  formulation  and 
evaluation  of  the  energetics  and  interatomic  Interactions  in  materials 
[7,37]  coupled  with  the  development  and  implementation  of  computational 
methods  and  simulation  techniques  [7,38],  open  new  avenues  for  investigations 
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of  the  microscopic  origins  of  complex  materials  phenomena.  In  particular 
large-scale  molecular  dynamics  computer  simulations,  which  are  in  a  sense 
computer  experiments,  where  the  evolution  of  a  system  of  interacting  particles 
is  simulated  with  high  spatial  and  temporal  resolution  via  direct  integration 
of  the  particles'  equations  of  motion,  have  greatly  enhanced  our  understanding 
of  a  broad  range  of  materials  phenomena. 

Although  our  knowledge  of  interfacial  processes  occurring  when  two 
material  bodies  are  brought  together  has  significantly  progressed  since  the 
original  presentation  by  Heinrich  Hertz  before  the  Berlin  Physical  Society  in 
January  1881  of  his  theory  of  the  contact  of  elastic  bodies  [17],  full  micro¬ 
scopic  understanding  of  these  processes  is  still  lacking.  Moreover,  it  has 
been  recognized  that  continuum  mechanics  is  not  fully  applicable  as  the  scale 
of  the  material  bodies  and  the  characteristic  dimension  of  the  contact  between 
them  are  reduced  [22,39].  Furthermore,  it  had  been  observed  [l8b,25]  that  the 
mechanical  properties  of  materials  exhibit  a  strong  dependence  on  the  size  of 
the  sample  (small  specimens  appear  to  be  stronger  than  larger  ones).  Since 
the  Junctions  between  contacting  solids  can  be  small,  their  mechanical 
properties  may  be  drastically  different  from  those  of  the  same  materials  in 
their  bulk  form.  Consequently,  the  application  of  the  newly  developed 
theoretical  and  experimental  techniques  to  these  problems  promises  to  provide 
significant  insights  concerning  the  microscopic  mechanisms  and  the  role  of 
surface  forces  in  the  formation  of  microcontacts  and  to  enhance  our 
understanding  of  fundamental  issues  pertaining  to  interfacial  adherence, 
microindentation,  structural  deformations,  and  the  transition  from  elastic  to 
elastoplastic  or  fully  developed  plastic  response  of  materials.  Additionally, 
studies  such  as  those  described  in  this  paper  allow  critical  assessment  of  the 
range  of  validity  of  continuum-based  theories  of  these  phenomena  and  could 
inspire  improved  analytical  formulations.  Finally,  knowledge  of  the 
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interactions  and  atomic-scale  processes  occurring  between  small  tips  and 
material  surfaces,  and  their  consequences,  is  of  crucial  importance  to 
optimize,  control,  interpret  and  design  experiments  employing  the  novel 
tip-based  microscopies  [6,7, 13.14, 16,  30,33-36,40-45]. 

To  this  end  we  investigate  in  this  paper  theoretically  [7]  and 
experimentally  [6]  the  interactions  between  a  metallic  tip  (Ni)  and  a  gold 
surface,  chosen  mainly  due  to  differences  in  their  mechanical  properties  such 
as  elastic  moduli,  yield,  hardness  and  strength  parameters  (e.g.,  the  elastic 
moduli  are  21  x  lO1^  N/m2  and  8.2  x  10^  N/m2  for  Ni  and  Au,  respectively 
[46]).  The  theoretical  studies  employed  molecular  dynamics  (MD)  simulations 
[7,16]  with  interatomic  interactions  described  via  the  many-body  potentials 
obtained  by  the  embedded-atom  method  (EAM)  [47]  which  have  been  recently 
applied  with  significant  success  in  studies  of  bulk  and  surface  properties  of 
a  number  of  metallic  systems  and  their  alloys.  The  experimental  measurements 
were  performed  using  AFM  configured  to  measure  the  force  between  a  tip  mounted 
on  a  cantilever  and  the  sample  surface  as  a  function  of  tip-to-sample 
separation  [6,45]. 

Our  theoretical  simulations  reveal  the  onset  of  an  instability  as  the  tip 
approaches  the  sample  causing  a  Jump-to-contact  (JC)  such  as  described  first 
[14]  via  calculations  employing  Lennard-Jones  (LJ)  potentials  and  further 
investigated  more  recently  using  different  potentials  for  nickel  [15]  and 
other  materials  [73.  He  find  that  for  our  system  the  JC  phenomenon  is 
associated  primarily  with  a  tip-induced  sample  deformation  which  begins  when 
the  distance  between  the  proximal  atomic  layers  of  the  two  interfacing 
materials  is  approximately  4.2  8  (i.e.,  at  a  separation  larger  than  the 
equilibrium  crystalline  interlayer  spacings),  and  that  the  process  involves 
large  atomic  displacements  (~2  8)  occurring  over  a  short  time  span  of  ~  1 
psec.  Furthermore,  we  discovered  that  lifting  the  tip  from  the  surface  after 
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contact  results  in  an  inelastic  deformation  of  the  sample  exhibiting  ductile 
extension  and  the  eventual  tear,  or  fracture,  of  the  topmost  Au  layer  which 
adheres  to  the  Ni  tip.  The  above  processes  are  portrayed  in  both  the  results 
of  the  simulations  and  measurements  as  a  marked  hysteresis  in  the  force  versus 
distance  relationship  recorded  along  the  axis  of  the  tip-sample  approach/ 
separation.  In  fact,  the  process  of  tear,  observed  in  the  simulations  during 
tip-sample  separation,  is  akin  to  mode-I  ductile  fracture  [31]  (i.e.,  load 
normal  to  the  fracture  plane).  Allowing  the  tip  to  further  advance  and 
penetrate  the  sample  surface  beyond  the  point  of  contact  indents  the  surface 
and  results  in  further  deformation  of  the  sample  characterized  by  an  adhesion- 
induced  flow  of  gold  atoms  which  wet  the  edges  of  the  Ni  tip,  generation  of 
slip  planes  in  the  Au  lattice  and  formation  of  point  defects.  Separating  the 
tip  and  sample  causes  the  sample  to  deform  ductilely  producing  an  extended 
crystalline  "neck"  that  stretches  between  the  sample  and  adherent  layers  on 
the  tip  until  the  neck  eventually  breaks.  Throughout  much  of  the  elongation 
process  the  neck  maintains  a  crystalline  layered  structure  while  reducing  in 
cross-sectional  area  as  it  extends.  The  elongation  mechanism  revealed  by  the 
simulations  consists  of  a  sequence  of  elastic  and  plastic  (yielding)  stages 
accompanied  by  atomic  structural  rearrangement.  Based  on  the  above  observa¬ 
tions  we  associate  the  calculated  and  measured  hystereses  in  the  force  versus 
distance  curves  with  the  formation,  stretching  and  breaking  of  bonds  due  to 
adhesion,  cohesion  and  decohesion,  and  with  inelastic  deformations  induced  by 
the  tlp-to-substrate  interaction. 

METHODOLOGY 

Prior  to  the  presentation  of  our  results  we  provide  pertinent  details  of 
our  studies,  noting  common  as  well  as  distinguishing  characteristics  between 
the  theoretical  and  experimental  modes  of  investigation. 
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A.  MD  Simulations 

MO  simulations  were  performed  for  an  Au  (001)  sample  containing  3  static 
(bottom)  and  8  dynamic  layers,  each  consisting  of  450  atoms  per  layer.  The 
sample  interacts  with  a  dynamic  Ni  tip  consisting  of  1400  atoms  arranged 
originally  as  a  pyramidal  (tapered)  tip  with  the  bottom  layer  (closest  to  the 
sample)  consisting  of  72  atoms  exposing  a  (001)  surface,  the  next  layer 
consisting  of  128  atoms  and  the  remaining  6  layers  containing  200  atoms  each. 
This  gives  the  tip  an  effective  radius  of  curvature  of  ~  30  X  (which  is 
approximately  50-100  times  smaller  than  the  tip  employed  in  the  experiment). 

In  addition  the  tip  interacts  with  a  static  Ni  holder  consisting  of  1176  atoms 
arranged  in  three  (001)  layers.  The  simulations  were  performed  at  300K  with 
temperature  control  imposed  only  on  the  deepest  dynamic  layer  of  the  Au  sample 
closest  to  the  static  layers.  No  significant  variations  in  temperature  were 
observed  during  the  simulations.  The  equations  of  motion  were  integrated 
using  a  5-th  order  predictor-corrector  algorithm  with  a  time  step  At  *  3.05 
fsec. 

As  aforementioned  the  interatomic  interactions  which  govern  the 
energetics  and  dynamics  of  our  system  are  modeled  via  the  embedded-atom  method 
(EAM)  [47]  which  has  been  applied  recently  with  significant  success  to  study 
equilibrium  and  nonequilibrium  properties  and  processes  in  metallic  elemental 
and  alloy  systems  [47-49],  In  this  method,  the  dominant  contribution  to  the 
cohesive  energy  of  the  material  is  viewed  as  the  energy  to  embed  an  atom  into 
the  local  electron  density  provided  by  other  atoms  of  the  system.  This 
background  density  is  determined  for  each  atom  as  the  superposition  of 
electronic  densities  from  the  other  atoms,  evaluated  at  the  location  of  the 
atom  in  question.  Thus,  the  total  cohesive  energy  is  represented  in  the  EAM 
by  a  many-body  embedding  functional,  supplemented  by  parametrized  short-range 
pair  interactions  due  to  inter-core  repulsion.  The  parameters  of  the  pair- 
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potentials  are  determined  via  fitting  to  a  number  of  bulk  equilibrium 
properties  of  the  metals  and  their  alloys,  such  as  lattice  constants,  heats  of 
sublimation,  elastic  constants,  vacancy-formation  energies  and  heats  of 
solution  [47]. 

Following  equilibration  of  the  system  at  300K  with  the  tip  outside  the 
range  of  interaction,  the  tip  was  lowered  slowly  toward  the  surface.  For  the 
initially  equilibrated  system  we  find  multilayer  relaxation  [50]  of  the 
Au(001)  surface,  whereby  the  first  ( topmost) -to-second  interlayer  spacing, 
d12,  contracts  by  7.5J  and  d^  expands  by  3.5%  relative  to  the  interlayer 
spacing  in  the  bulk.  The  layer  relaxation  at  the  surface  of  the  Ni  tip  is 
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insignificant  and  at  equilibrium  the  sides  of  the  tapered  part  of  the  tip 

expose  small  (111)  facets.  The  calculated  surface  energies  (at  OK)  for  Ni 
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(001)  and  Au  (001)  are  1657  mJ/m  and  964  mj/m  ,  respectively,  in  close 
agreement  with  calculations  [47a]  employing  a  slightly  different 
parameterization  of  the  EAM  potentials. 

Motion  of  the  tip  occurs  by  changing  the  position  of  the  tip-holder 
assembly  in  increments  of  0.25  S  over  500  At.  After  each  increment  the  system 
is  fully  relaxed,  i.e.,  dynamically  evolved,  until  no  discernable  variations 
in  system  properties  are  observed  beyond  natural  fluctuations. 

B.  AFM  Measurements 

Forces  may  be  measured  experimentally  using  AFM  [6,331;  both  attractive 
and  repulsive  forces  can  be  measured  as  well  as  the  adhesive  force  necessary 
to  separate  the  tip  and  sample  surface  after  contact.  Fig.  1  depicts 
schematically  the  forces  acting  between  the  tip  (which  is  mounted  on  a 
cantilever  beam)  and  sample  as  a  function  of  the  separation  D  between  the 
cantilever  tip  and  sample.  The  arrows  are  used  to  guide  the  eye  throughout 
the  full  interaction  cycle.  The  cycle  starts  with  the  sample  far  away  and  the 
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cantilever  in  its  rest  position.  As  D  decreases,  the  cantilever  bends  towards 
the  sample  such  that  at  any  equilibrium  separation  D  the  attractive  force,  F, 
balances  the  restoring  force  of  the  cantilever  defined  by  its  effective  spring 
constant  k  times  its  deflection.  However,  if  the  magnitude  of  the  gradient 
cf  the  attractive  force  dF/dD  exceeds  k  (point  A),  the  cantilever  will  Jump 
into  contact  with  the  sample  (point  A').  This  instability  is  governed  by  the 
stiffness  of  the  cantilever  beam  relative  to  the  long-range  forces,  while  in 
the  HD  simulations,  where  the  cantilever  is  modeled  by  a  rigid  tip-holder 
(i.e.,  an  infinitely  stiff  cantilever  beam),  the  Jump- to -con tact  instability 
is  driven  by  the  inherent  stiffness  (related  to  the  cohesive  strength)  of  the 
tip  and  substrate  materials. 

On  reversing  the  direction  of  the  sample,  the  cantilever  will  Jump  away 
from  the  sample  at  B  to  some  point  B'  giving  rise  to  hysteresis  in  the 
measured  force  curve,  the  magnitude  of  which  depends  on  k  and  F.  Thus  in  the 
experiment,  the  degree  of  resolution,  i.e.,  the  ability  to  track  the  force 
versus  distance  curve  for  all  tip- to-sur face  separations,  depends  on  the 
selection  of  the  cantilever.  For  example,  if  k  is  at  any  time  less  than 
dF/dD,  the  dotted  segment  of  the  interaction  force  curve  Dft  <  D  <  D0  (see 
Fig.  1)  is  inaccessible.  However,  if  k  is  always  greater  than  dF/dD,  as  for 
our  experiment,  the  cantilever  dependent  instability  can  be  practically 
eliminated,  thus  enabling  a  faithful  measurement  of  the  consequences  of  the 
tip-to-sample  Interatomic  interactions  and  the  resultant  hysteresis,  due  to 
adhesion,  in  the  force  versus  separation  curve. 

The  AFM  instrument  employs  a  cantilever  beam  made  from  Ni  wire  (0.25  mn 
dia.)  and  bent  into  the  shape  of  a  'L'  whose  long  and  short  dimensions  are 
approximately  6  and  2  no,  respectively.  The  short  part  of  the  bent  wire  was 
chemically  etched  with  hydrochloric  acid  into  a  tip  with  a  radius  of  curvature 
of  ~  200  nm  as  determined  by  scanning  electron  microscopy.  The  cantilever's 
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spring  constant  is  calculated  to  be  ~  5000  N/m.  Its  deflection  is  measured 
with  a  tunneling  microscope.  Details  of  the  instrumentation  are  described 
elsewhere  [6].  Whereas  in  the  simulation  the  tip  is  moved,  in  the  experiment 
the  sample  is  moved  by  means  of  a  piezoelectric  actuator  at  a  rate  of  5 
nm/sec.  Sample  velocity  was  chosen  to  be  intermediate  to  the  thermal  drift 
rate  (<  1  nm/min)  and  tunneling  microscope  slew  rate  (~  100  nm/sec). 

The  AFM  measurements  were  done  in  a  dry  box  under  dry  nitrogen  using  tips 
and  samples  that  were  prepared  in  air  and  then  quickly  transferred  to  the  dry 
box.  The  sample  was  an  evaporated  Au  film  (approximately  100  nm  thick)  that 
was  cleaned  in  a  solution  of  sodium  dichromate  and  concentrated  sulfuric  acid 
and  rinsed  with  distilled  water  until  the  surface  was  hydrophilic. 

RESULTS  AND  DISCUSSION 

Measured  and  simulated  force  versus  distance  curves  are  shown  in  Figs.  2 
and  3,  respectively,  as  well  as  the  calculated  potential  energy  versus 
distance  (Fig.  3c).  In  both  cases  results  for  tip-to-sample  approach  followed 
by  separation  are  shown,  for  adhesive  contact  (Figs.  2a  and  3a)  and 
indentation  (Figs.  2b  and  3b, 3c)  studies.  As  discussed  above,  the  simulations 
correspond  to  a  case  of  a  rigid  cantilever  and  therefore  the  recorded 
properties  of  the  system  as  the  tip-holder  assembly  approaches  or  retracts 
from  the  sample  portray  directly  consequences  of  the  interatomic  interactions 
between  the  tip  and  the  sample.  The  distance  scale  that  we  have  chosen  in 
presenting  the  calculated  results  is  the  separation  (denoted  as  dha)  between 
the  rigid  (static)  holder  of  the  tip  and  the  static  gold  lattice  underlying 
the  dynamic  substrate.  The  origin  of  the  distance  scale  is  chosen  such  that 
d^a  s  0  after  Jump-to-contact  occurs  (d^a  >  0  when  the  system  is  not  advanced 
beyond  the  JC  point  and  dha  <  0  corresponds  to  indentation).  Since  the 
dynamic  Ni  tip  and  Au  substrate  atoms  displace  in  response  to  the  interaction 
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between  them,  the  distance  d^  does  not  give  directly  the  actual  separation 
between  regions  in  the  dynamic  tip  and  substrate  material.  The  actual 
relative  distances.  dfca,  between  the  bottom  part  of  the  tip  (averaged 
z-position  of  atoms  in  the  bottom-most  layer  of  the  tip)  and  the  surface 
(averaged  z-position  of  the  topmost  layer  of  the  Au  surface,  calculated  for 
atoms  in  the  first  layer  away  from  the  perturbed  region  in  the  vicinity  of  the 
tip)  are  given  by  the  letter  symbols  in  Fig.  3(a,b)  Note  that  the  distance 
between  the  bottom  of  the  tip  and  the  gold  atoms  in  the  region  inmed lately 
underneath  it  may  differ  from  dfca.  Thus  for  example  when  dhg  =  0  (point  D  in 
Figs.  3a, b)  the  tip  to  unperturbed  gold  distance,  dfca,  is  3.8  8,  while  the 
average  distance  between  the  bottom  layer  of  the  tip  and  the  adherent  gold 
layer  in  inmed late  contact  with  it  is  2.1  8. 

A.  Tip-Sample  Approach 

Comparison  of  Figs.  2a  and  3a  reveals  similarity  between  the  measured  and 
calculated  curves  showing  a  monotonic  increase  in  the  magnitude  of  the 
attractive  force  as  the  tip  approaches  the  sample  and  a  hysteresis  during 
separation.  Note,  however,  that  in  the  experiment  (Fig.  2a)  the  magnitude  of 
the  force  and  the  distance  from  the  surface  at  which  it  begins  to  deviate  from 
zero  are  much  larger  than  in  the  calculated  data  (Fig.  3a).  These  differences 
are  due  primarily  to  differences  in  tip  size  and  to  the  neglect  of  long-range 
interactions  (e.g.,  van  aer  Waals  forces)  in  the  calculations.  Also,  in  the 
experiment,  it  is  difficult  to  determine  where  the  cantilever  tip  is  in 
relation  to  the  surface  until  contact  has  already  been  made.  Therefore  it  is 
likely  that  a  slight  indentation  has  occurred  in  the  data  shown  in  Fig.  2a. 
This  problem  can  be  alleviated  in  future  experiments  by  using  phase  sensitive 
detection  to  denote  contact  and  reverse  the  direction  of  the  sample  motion. 
Other  differences  between  theory  and  experiment  include  tip  or  sample 
roughness  and  exposure  to  air  during  sample  preparation.  We  stress  that  in 
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our  experiments  cantilever-dependent  instability  has  been  ruled  out  as  a 
possible  origin  of  the  observed  hysteresis.  Furthermore,  hysteresis  in  the 
piezoelectric  actuators  used  in  the  instrument  is  negligible.  Rather,  the 
hysteresis  originates  from  dynamical  processes  induced  by  the  tip-substrate 
adhesive  interactions  as  revealed  by  the  simulations. 

Following  an  initial  slow  variation  of  the  force  between  the  Au  substrate 
and  the  Ni  tip  we  observe  in  the  simulations  the  onset  of  an  instability, 
signified  by  a  sharp  increase  in  the  attraction  between  the  two  (see  Fig.  3a 
as  well  as  Figs.  3b  and  3c  where  the  segments  corresponding  to  lowering  of  the 
tip  up  to  the  point  D  describe  the  same  stage  as  that  shown  in  segment  AD  in 
Fig.  3a.)  which  is  accompanied  by  a  marked  decrease  in  the  potential  energy  of 
the  system  (see  sudden  drop  of  Ep  in  Fig.  3c  as  dhg  approaches  zero  from  the 
right).  We  note  the  rather  sudden  onset  of  the  instability  which  occurs  only 
for  separations  smaller  than  0.25  R  (marked  by  an  arrow  on  the  curve  in 
Fig.  3a).  Our  simulations  reveal  that  in  response  to  the  imbalance  between 
the  forces  on  atoms  in  each  of  the  materials  and  those  due  to  intermetal lie 
interactions  a  Jump-to-contact  (JC)  phenomenon  occurs  via  a  fast  process  where 
Au  atoms  in  the  region  of  the  surface  under  the  Ni  tip  displace  by 
approximately  2  R  toward  the  tip  in  a  short  time  span  of  ~  1  ps  (see  bulging 
of  the  gold  surface  shown  in  Fig.  4a  where  the  atomic  configuration  after  the 
JC  is  depicted).  After  the  Jump-to-contact  occurs  the  distance  between  the 
bottom  layer  of  the  Ni  tip  and  the  layer  of  adherent  Au  atoms  in  the  region 
immediately  underneath  it  decreases  to  2.1  R  from  a  value  of  4.2  R.  In 
addition  to  the  adhesive  contact  formation  between  the  two  surfaces  an 
adhesion-induced  partial  wetting  of  the  edges  of  the  Ni  tip  by  Au  atoms  is 
observed  (see  Fig.  4a). 

The  Jump-to-contact  phenomenon  in  metallic  systems  is  driven  by  the 
marked  tendency  of  the  atoms  at  the  interfacial  regions  of  the  tip  and 
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substrate  materials  to  optimise  their  embedding  energies  (which  are  density 
dependent,  deriving  from  the  tails  of  the  atomic  electronic  charge  densities) 
while  maintaining  their  individual  material  cohesive  binding  (in  the  Ni  and 
Au)  albeit  strained  due  to  the  deformation  caused  by  the  atomic  displacements 
during  the  JC  process. 

Further  insight  into  the  JC  process  is  provided  by  the  local  hydrostatic 

pressure  in  the  materials  (evaluated  as  the  trace  of  the  atomic  stress  tensors 

[51])  shown  in  Fig.  5a  after  contact  formation  (i.e.,  point  D  in  Fig.  3a). 

The  pressure  contours  reveal  that  atoms  at  the  periphery  of  the  contact  zone 

e 

(at  X  s  ±0.19  and  Z  =  0.27)  are  under  extreme  tensile  stress  (-  10  atm  = 

10  2 

-  10  N/m  *  -lOGPa).  In  fact  we  observe  that  the  tip  as  well  as  an  extended 
region  of  the  substrate  in  the  vicinity  of  the  contact  zone  are  under  tension. 
Both  the  structural  deformation  profile  of  the  system  and  the  pressure 
distribution  which  we  find  in  our  atomistic  HD  simulations  are  similar,  in 
general  terms,  to  those  described  by  certain  modern  contact  mechanics  theories 
[19-22]  where  the  influence  of  adhesive  interactions  is  Included. 

B.  Tip-Substrate  Separation  After  Contact 

Starting  from  contact  the  force  versus  distance  (Fz  vs.  d^)  curve 
exhibits  a  marked  hysteresis  seen  both  experimentally  (Fig.  2a)  and 
theoretically  (Fig.  3a)  as  the  surfaces  are  separated.  We  remark  that,  in 
the  simulation  and  the  measurements,  separating  the  surfaces  prior  to  contact 
results  in  no  hysteresis.  The  hysteresis  is  a  consequence  of  the  adhesive 
bonding  between  the  two  materials  and,  as  demonstrated  by  the  simulation, 
separation  is  accompanied  by  inelastic  processes  in  which  the  topmost  layer  of 
the  Au  sample  adheres  to  the  Ni  tip.  (See  the  configuration  shown  in  Fig.  4b 
which  corresponds  to  the  distance  dfcg  s  7.5  8  marked  J  in  Fig.  3a).  The 
mechanism  of  the  process  is  demonstrated  by  the  pressure  contours  during 
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liftoff  of  the  tip  shown  in  Fig.  5b,  recorded  for  the  configuration  marked  G 
(d.  =  5.5  8  in  Fig.  3a).  As  seen  the  maximum  tensile  stress  is  located  near 
the  edges  of  the  adhesive  contact.  We  further  observe  that  the  diameter  of 
the  contact  area  decreases  during  lifting  of  the  tip,  resulting  in  the 
formation  of  a  thin  "adhesive  neck"  due  to  ductile  extension,  which  stretches 
as  the  process  continues,  ultimately  breaking  at  a  distance  dfc3  of  ~  9-10  8. 
The  evolution  of  adhesion  and  tear  mechanisms  which  we  observe  can  be 
classified  as  mode-I  fracture  [31],  reemphasizing  the  importance  of  forces 
operating  across  the  crack  in  modeling  crack  propagation  [31,32]. 

C.  Indentation 

We  turn  now  to  theoretical  and  experimental  results  recorded  when  the  tip 
is  allowed  to  advance  past  the  Jump- to-con tact  point,  i.e.,  indentation  (see 
Figs.  2b,  3(b,c),  and  Me).  As  evident  from  Fig.  (3b),  decreasing  the 
separation  between  the  tip  and  the  substrate  causes  first  a  decrease  in  the 
magnitude  of  the  force  on  the  tip  (i.e.,  less  attraction,  see  segment  DL)  and 
an  increase  in  the  binding  energy  (i.e.,  larger  magnitude  of  the  potential 
energy,  shown  in  Fig.  3c).  However,  upon  reaching  the  point  marked  L  in  Fig. 
(3b)  a  sharp  increase  in  the  attraction  occurs,  followed  by  a  monotonic 
decrease  in  the  magnitude  of  the  force  till  F  =  0  (point  M  in  Fig.  3b)  at  d,. 

=  0.8  8.  The  variations  of  the  force  (in  the  segment  DLM)  are  correlated  with 
large  deformations  of  the  Au  substrate  (see  the  atomic  configuration  in  Fig. 
4c,  corresponding  to  point  M  in  Fig.  3b).  In  particular,  the  nonmonotonic 
feature  (near  point  L)  results  from  tip-induced  flow  of  gold  atoms  which 
relieve  the  increasing  stress  via  wetting  of  the  sides  of  the  tip.  Indeed  the 
atomic  configurations  (Figs,  4c  and  5c)  display  a  "piling-up"  around  the  edges 
of  the  indenter  due  to  atomic  flow  driven  by  the  deformation  of  the  Au 
substrate  and  the  adhesive  interactions  between  the  Au  and  Ni  atoms.  Further 
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indentation  is  accompanied  by  slip  of  Au  layers  (along  (111)  planes)  and  the 
generation  of  interstitial  defects  (see  atomic  trajectories  and  atomic 
configuration,  corresponding  to  point  M,  in  Figs.  4c  and  5c,  respectively). 

In  addition,  the  calculations  predict  that  during  the  indentation  process  a 
small  number  of  Ni  atoms  diffuse  into  the  surrounding  Au,  occupying 
substitutional  sites.  Furthermore  the  calculated  pressure  contours  at  this 
stage  of  indentation,  shown  in  Fig.  5d,  demonstrate  that  the  substrate  surface 
zone  in  the  vicinity  of  the  edges  of  the  tip  is  under  tensile  stress,  while 
the  deformed  region  under  the  tip  is  compressed  with  the  maximum  pressure 
(8.2GPa)  occurring  at  about  the  fifth  Au  layer  below  the  center  of  the  Ni 
tip-indenter .  The  general  characteristics  of  the  pressure  (and  stress) 
distributions  obtained  in  our  indentation  simulations  correspond  to  those 
associated  [11,19,23]  with  the  onset  and  development  of  plastic  deformation  in 
the  substrate. 

Experimentally,  advancing  the  sample  past  the  contact  point  is  noted  by 
the  change  in  slope  of  the  force  as  the  increasing  repulsive  forces  push  the 
tip  and  cantilever  back  towards  their  rest  position,  as  shown  in  Fig.  2b.  Ue 
remark  that  the  calculated  pressures  from  the  simulations  compare  favorably 
with  the  average  contact  pressure  of  ~  3  GPa  determined  experimentally  [45]  by 
dividing  the  measured  attractive  force  by  the  estimated  circular  contact  area 
of  radius  20  nm. 

D.  Tip-Substrate  Separation  After  Indentation 

Reversal  of  the  direction  of  the  tip  motion  relative  to  the  substrate 
from  the  point  of  zero  force  (point  M  in  Fig.  3b)  results  in  the  force-  and 
potential  energy-  versus  distance  curves  shown  in  Figs.  3b  and  3c.  The  force 
curve  exhibits  first  a  sharp  monotonic  increase  in  the  magnitude  of  the 
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attractive  force  (segment  MN  in  Fig.  3b)  with  a  corresponding  increase  in  the 
potential  energy  (Fig.  3c).  During  this  stage  the  response  of  the  system  is 
mostly  elastic  accompanied  by  the  generation  of  a  small  number  of  vacancies 
and  substitutional  defects  in  the  substrate.  Past  this  stage  the  force  and 
energy  curves  versus  tip-to-sample  separation  exhibit  a  nonmonotonic  behavior 
which  is  associated  mainly  with  the  process  of  elongation  of  the  connective 
neck  which  forms  between  the  substrate  and  the  retracting  tip. 

To  illustrate  the  neck  formation  and  elongation  process  we  show  in  Fig.  6 
a  sequence  of  atomic  configurations  corresponding  to  the  maxima  in  the  force 
curve  (Fig.  3b,  points  marked  0,  Q,  S,  U,  W  and  X).  As  evident,  upon 
increased  separation  between  the  tip-holder  and  the  substrate  a  connective 
neck  forms  consisting  mainly  of  gold  atoms  (see  atomic  configurations  shown  in 
Fig.  4d  and  4e).  The  mechanism  of  elongation  of  the  neck  involves  atomic 
structural  transformations  whereby  in  each  elongation  stage  atoms  in  adjacent 
layers  in  the  neck  disorder  and  then  rearrange  to  form  an  added  layer,  i.e.,  a 
more  extended  neck  of  a  smaller  cross-sectional  area.  Throughout  the  process 
the  neck  maintains  a  layered  crystalline  structure  (see  Figs.  6  and  4e)  except 
for  the  rather  short  structural  transformation  periods,  corresponding  to  the 
sharp  variations  in  the  force  curve,  (see  segments  PQ,  RS,  TU  and  VW  in  Fig. 
3b)  and  the  associated  features  in  the  calculated  potential  energy  shown  in 
Fig.  3c  where  the  minima  correspond  to  ordered  layered  structures  after  the 
structural  rearrangements.  We  note  that  beyond  the  initial  formation  stage, 
the  number  of  atoms  in  the  connective  neck  region  remains  roughly  constant 
throughout  the  elongation  process. 

Further  insight  into  the  microscopic  mechanism  of  elongation  of  the 
connective  neck  can  be  gained  via  consideration  of  the  variation  of  the  second 
invariant  of  the  stress  deviator,  which  is  related  to  the  von  Mises  shear 
strain-energy  criterion  for  the  onset  of  plastic  yielding  [ 18a ,19].  Returning 
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to  the  force  and  potential  energy  curves  shown  in  Figs.  3b  and  :>c,  we  have 
observed  that  between  each  of  the  elongation  events  (i.e.t  layer  additions, 
points  marked  Q,  S,  U,  W  and  X)  the  initial  response  of  the  system  to  the 
strain  induced  by  the  increased  separation  between  the  tip-holder  and  the 
substrate  is  mainly  elastic  (segments  OP,  QR,  ST,  UV  in  big.  3b,  and 
correspondingly  the  variations  in  Fig.  3c),  accompanied  by  a  gradual  increase 
of  '/J. £ ,  and  thus  the  stored  strain  energy.  The  onsets  of  the  stages  of 
structural  rearrangements  are  found  to  be  correlated  with  a  critical  maximum 
value  of  vCT,  of  about  3  GPa  ( occur ing  for  states  at  the  end  of  the  intervals 
marked  OP,  QR,  ST  and  UV  in  Fig.  3b)  localized  in  the  neck  in  the  region  of 
the  ensuing  structural  transformation.  After  each  of  the  elongation  events 
the  maximum  value  of  VT,  (for  the  states  marked  Q,  S,  U,  W  and  X  in  Fig.  3b) 
drops  to  approximately  2  GPa. 

In  this  context,  it  is  interesting  to  remark  that  the  value  of  the  normal 
component  of  the  force  per  unit  area  in  the  narrowest  region  of  the  neck 
remains  roughly  constant  (~  lOGPa)  throughout  the  elongation  process, 
increasing  by  about  20£  prior  to  each  of  the  aforementioned  structural 
rearrangements.  This  value  has  been  estimated  both  by  using  the  data  given  in 
Figs.  3c  and  the  cross  sectional  areas  from  atomic  configuration  plots  (such 
as  given  in  Fig.  6),  and  via  a  calculation  of  the  average  axial  component  (zz 
element)  of  the  atomic  stress  tensors  [51]  in  the  narrow  region  of  the  neck. 

We  note  that  the  above  observations  constitute  atomic-scale  realizations  of 
basic  concepts  which  underlie  macroscopic  theories  of  materials  behavior  under 
load  [17-19]. 

A  typical  distribution  of  the  stress,  prior  to  a  structural  trans¬ 
formation  is  shown  in  Fig.  7  (shown  for  the  state  corresponding  to  the  point 
marked  T  in  Fig.  3b).  As  seen,  the  maximum  of  is  localized  about  a  narrow 
region  around  the  periphery  in  the  strained  neck.  Comparison  between  the 
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atomic  configuration  at  this  stage  (see  Fig.  7,  or  the  very  similar  configura¬ 
tion  shown  in  Fig.  6c)  and  the  configuration  after  the  structural  trans¬ 
formation  has  occurred  (see  Fig.  6d,  corresponding  to  the  point  marked  U  in 
Fig.  3b)  illustrates  the  elongation  of  the  neck  by  the  addition  of  a  layer  and 
accompanying  reduction  in  areal  cross  section.  We  note  that  as  the  height  of 
the  connective  neck  increases  the  magnitude  of  the  variations  in  the  force  and 
potential  energy  during  the  elongation  stages  diminishes.  The  behavior  of  the 
system  past  the  state  shown  in  Fig.  6f  (corresponding  to  the  point  marked  X  in 
the  force  curve  shown  in  Fig.  3b)  is  similar  to  that  observed  at  the  final 
stages  of  separation  after  jump-to-contact  (Fig.  3b) ,  characterized  by  strain 
induced  disordering  and  thinning  in  a  narrow  region  of  the  neck  near  the  gold 
covered  bottom  of  the  tip  and  eventual  fracture  of  the  neck  (occuring  for  a 
tip-to-substrate  distance  dfc3  a  18  8),  resulting  in  a  Ni  tip  whose  bottom  is 
covered  by  an  adherent  Au  layer. 

The  theoretically  predicted  increased  hysteresis  upon  tip-substrate 
separation  following  indentation,  relative  to  that  found  after  contact 
(compare  Figs.  3a  and  3b),  is  also  observed  experimentally,  as  shown  in  Fig. 
2b.  In  both  theory  and  experiment  the  maximum  attractive  force  after 
indentation  is  roughly  50%  greater  than  when  contact  is  first  made.  Note 
however  that  the  nonmonotonic  features  found  in  the  simulations  (Fig.  3b)  are 
not  discernible  in  the  experiment  which  is  apparently  not  sufficiently 
sensitive  to  resolve  such  individual  atomic-scale  events  when  averaged  over 
the  entire  contact  area. 

Finally  we  remark  that  scanning  Auger  microscopy  (SAM)  was  used  to  detect 
gold  transfer  onto  the  nickel  tip  following  indentation.  The  lack  of 
sensitivity  of  SAM  to  concentrations  of  materials  less  than  0. If  prohibits  us 
from  detecting  Au  transfer  onto  the  tip  for  the  typical  small  loads  used  in 
the  measurements  discussed  in  this  paper.  However,  we  did  observe  the 
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presence  of  Au  on  the  Nl  tip  following  the  application  of  a  load  of  ~  100  mN 
between  a  freshly  cleaned  Ni  tip  and  a  Au  surface,  just  moments  before  placing 
the  tip  into  the  Auger  vacuum  chamber  for  analysis. 


CONCLUDING  REMARKS 

The  recent  emergence  and  proliferation  of  proximal  probes,  in  particular 
tip-based  microscopies  [33-36],  sensitive  to  nano-  and  subnano-  meter  scale 
structures  provides  compelling  opportunities  for  studies  of  these  structures 
which  are  key  to  the  science  base  of  many  venerable  technological  problems 
[36].  In  addition  these  probes,  coupled  with  advances  in  theoretical 
understanding  of  the  energetics  and  interaction  mechanisms  in  materials  and 
the  development  of  computer-based  materials  modeling  and  simulation 
techniques,  open  new  avenues  for  exploration  of  new  scientific  concepts  and 
novel  materials  properties  and  processes  on  the  subnanometer  scale  [36,52], 

In  this  paper  we  have  presented  results  of  Joint  theoretical,  molecular 
dynamics  simulations,  and  experimental,  atomic  force  microscopy,  studies  of 
the  mechanisms  and  properties  of  interroetallic  adhesive  interactions,  contact 
formation,  nanoindentation  and  mechanical  response.  Our  studies  show  that 
contact  formation  between  a  hard  tip  {nickel)  approaching  a  soft  metallic 
substrate  (gold)  is  associated  with  an  atomic-scale  instability  which  leads  to 
a  Jump-to-contact  phenomenon  which  involves  an  inelastic  response  of  the  atoms 
in  the  proximal  interfacial  region  of  the  gold  substrate.  Indentation  of  the 
surface  by  advancing  the  tip  beyond  the  point  of  contact  results  in  the  onset 
of  plastic  yielding,  adhesion- induced  atomic  flow  and  generation  of  slip  in 
the  surface  region  of  the  gold  substrate.  Separating  the  two  materials  from 
contact  results  in  adhesion- induced  wetting  of  the  tip  by  gold  atoms, 
adherence  of  a  gold  monolayer  to  the  nickel  tip,  formation  of  an  atomically 
thin  connective  neck  and  eventual  fracture.  Furthermore,  retracting  the  tip 
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from  the  sample  after  indentation  results  in  ductile  extension  of  the 
substrate  and  formation  of  a  connective  crystalline  neck  which  elongates, 
while  reducing  in  cross-sectional  area,  via  structural  transformations 
involving  elastic  and  yielding  stages. 

The  above  microscopic  processes  are  portrayed  in  the  theoretically 
calculated  and  experimentally  measured  force  versus  distance  curves  (Figs.  2 
and  3),  which  exhibit  pronounced  hysteresis  upon  tip-to-sample  approach  and 
subsequent  separation.  We  note  however,  that  due  to  experimental  constraints 
(such  as  the  need  to  employ  a  cantilever  of  finite  flexibility,  sample  and  tip 
cleanliness,  and  tip  size)  certain  characteristic  features  predicted  by  the 
simulations  (i.e.,  the  nonmonotonic  behavior  of  the  force  vs.  distance 
relationship  during  separation  after  indentation,  associated  with  atomic 
structural  transformations  during  elongation  of  the  connective  neck)  could  not 
be  resolved  in  the  present  measurements. 

Our  Investigations  provide  the  impetus  for  further  combined  theoretical 
and  experimental  investigations  of  the  microscopic  mechanisms  of  adhesion, 
contact  formation  and  atomic-scale  mechanical  response  processes  in  materials, 
motivating  a  critical  assessment  of  the  range  of  validity  of  continuum 
theories  and  reformulation  of  contact  mechanics  formalisms  [10,19,531  to 
incorporate  an  atomistic  description  of  the  processes  of  adhesion,  deforma¬ 
tion,  wetting  and  fracture.  In  addition,  the  results  presented  here  are 
pertinent  to  the  general  issue  of  the  consequences  of  tip-substrate  inter¬ 
actions  in  tip  based  microscopies  (STM  and  AFM)  and  for  studies  of  the 
transition  from  tunneling  to  point  contact  in  STM  [12,13,54-56].  Finally, 
our  studies  suggest  a  method  for  preparation  of  atomic  size  contacts  and 
atomically  thin  wires  (via  the  process  of  contact  formation  between  a  metallic 
tip  and  a  soft  metal  substrate,  followed  by  gentle  separation)  which  could  be 
used  to  study  quantum  conductance  in  narrow  constrictions  [57]. 
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Figures 

Fig.  1.  (a)  Schematic  diagram  showing  the  relative  positions  of  the  tunneling 

tip,  cantilever  tip  and  sample.  The  tunneling  tip  is  under  feedback 
control,  operating  in  the  constant  current  mode.  The  piezos  are  used 
to  move  the  tunneling  tip  and  sample.  The  cantilever  for  this 
experiment  was  made  from  a  piece  of  Ni  wire  bend  into  the  3hape  of  a 
'L'  with  a  spring  constant  of  5000  N/m.  The  short  part  of  the  ' L ' 
was  chemically  etched  into  a  tip  with  a  radius  of  curvature  of 
approximately  200  nm.  (b)  Schematic  diagram  illustrating  the  force 
measurement  technique. 

Fig.  2.  Experimentally  measured  force  versus  tip-to-sample  distance  relation¬ 
ship  between  a  Ni  tip  and  Au  sample  for  (a)  contact  followed  by 
separation  and  (b)  indentation  followed  by  separation.  The  force 
curves  were  obtained  using  AFM  in  dry  nitrogen. 

Fig.  3.  Calculated  force,  Fz,  versus  tip-to-sample  distance,  d^,  relation¬ 
ship  between  a  Ni  tip  and  an  Au  sample  for:  (a)  approach  and  Jump- 
to-contact  followed  by  separation;  (b)  approach,  Jump-to-contact, 
indentation  and  subsequent  separation,  d^  denotes  the  distance 
between  the  rigid  tip-holder  assembly  and  the  static  substrate  of  the 
Au  surface  (d^  =  0  at  the  Jump-to-contact  point,  marked  D) .  The 
capital  letters  on  the  curves  denote  the  actual  distances,  dfcg, 
between  the  bottom  part  of  the  Ni  tip  and  the  Au  surface;  in  (a): 

A  =  5.7  8,  B  =  5.2  8,  C  =  4.7  8,  D  =  3.8  8,  E  =  4.4  8,  F  =  4.85  8, 

G  =  5.5  8,  H  =  5.9  8,  I  =  6.2  8  and  J  =  7.5  8;  in  (b):  D  =  3.8  8, 

L  =  2.4  8,  M  =  0.8  8,  N  =  2.6  8  0  =  3-0  8,  P  =  3.8  8,  Q  =  5.4  8, 

R  =  6.4  8,  S  =  7.0  8,  T  =  7.7  8,  U  =  9.1  8,  V  =  9.6  8,  W  =  10.5  8  and 

X  s  12.8  8.  (c)  Potential  energy  of  the  system  for  a  complete  cycle 

of  the  tip  approach,  Jump-to-contact,  indentation  and  subsequent 
separation.  Forces  in  units  of  nN,  energy  in  eV  and  distances  in  8. 
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Fig.  4.  Atomic  configurations  generated  by  the  MD  simulations,  (a)  after 
jump-to-contact  (see  point  D  in  Fig.  3a).  Note  bulging  of  the  Au 
substrate  under  the  Ni  tip  and  partial  wetting  of  the  tip  edges. 

(b)  Separation  after  contact  (point  J  in  Fig.  3a)  illustrating  ad¬ 
herence  of  the  top  Au  layer  to  the  Ni  tip  and  the  formation  of  an 
atomically  thin  connective  neck,  (c)  A  cut  through  the  system  at 
the  point  of  maximum  indentation  (M  in  Fig.  3b),  illustrating  de¬ 
formation  of  the  Au  substrate  and  slip  along  (111)  planes  of  the  sub¬ 
strate,  (d)  separation  after  indentation  (point  X  in  Fig.  3b), 

illustrating  wetting  of  the  Ni  tip  by  Au  atoms,  faceting  of  the  Ni 

( 

tip  exposing  (111)  facet  planes,  incorporation  of  Ni  atoms  in 
substitutional  sites,  and  formation  of  an  extensive  connective  neck 
between  the  tip  and  the  substrate.  Note  the  crystalline  character  of 
the  neck  and  the  incorporation  of  atoms  from  the  first,  second  and 
third  topmost  layers  of  the  Au  surface  as  well  as  several  Ni  atoms, 
(e)  A  cut  through  the  system  shown  in  (d)  illustrating  the  crystal¬ 
line  structure  of  the  neck  and  the  extent  of  structural  deformations 
of  the  substrate.  In  these  figures  the  Ni  tip  atoms  are  colored  red; 
atoms  of  the  top  layer  of  the  Au  surface  in  yellow,  second  layer  in 
blue,  third  layer  in  green,  fourth  in  yellow,  etc. 

Fig.  5.  Calculated  pressure  contours  and  atomic  configurations  viewed  along 
the  [010]  direction,  in  slices  through  the  system.  The  Ni  tip 
occupies  the  topmost  eight  atomic  layers.  Short-time  atomic  trajec¬ 
tories  appear  as  dots.  Distance  along  the  X  and  Z  directions  in 
units  of  X  =  1  and  Z  =  1  corresponding  to  61.2  ft  each.  Solid 
contours  correspond  to  tensile  stress  (i.e.,  negative  pressure)  and 
dotted  ones  to  compressive  stress,  (a)  after  Jump-to-contact  (point 
D  in  Fig.  3a,  see  also  Fig.  4a).  The  maximum  magnitude  of  the 
tensile  (i.e.,  negative  pressure,  lOGPa,  is  at  the  periphery  of  the 
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contact,  (X,Z)  =  (±  0.19,0.27).  The  contours  are  spaced  with  an 

increment,  A,  of  1  GPa.  Thus  the  contours  marked  e,  f  and  g 

correspond  to  -6,  -5  and  -4  GPa,  respectively,  (b)  During  separation 

following  contact,  (point  G  in  Pig.  3a).  The  maximum  tensile 

pressure  (marked  a),  — 9GPa,  is  at  the  periphery  of  tr.e  contact  at 

(X,Z)  equal  to  (0.1,0.25)  and  (-0.04,0.25).  A  =  0.9  GPa.  The  marked 

contours  h,  i,  j  and  k  correspond  to  -2.5,  -1.6,  -0.66  and  0.27GPa, 

respectively,  (c)  Short-time  particle  trajectories  at  the  final 

stage  of  relaxation  of  the  system,  corresponding  to  point  H  in  Fig. 

3b,  (i.e.,  F  =0).  Note  slip  along  the  [111]  planes  in  the 
z 

substrate,  (d)  Pressure  contours  corresponding  to  the  final 
configuration  shown  in  (c).  Note  the  development  of  compressive 
pressure  in  the  substrate  which  maximizes  in  the  region  of  the 
contour  marked  l  (8.2GPa).  The  increment  between  contours  A  = 

1.4GPa.  The  contours  marked  a  and  e  correspond  to  -6.4GPa  and 
-I.IGPa,  respectively,  and  those  marked  f  and  g  to  0.2  and  1.6GPa. 

Fig.  6.  Atomic  configurations  in  slices  through  the  system  illustrating 
the  formation  of  a  connective  neck  between  the  Ni  tip  and  the  Au 
substrate  during  separation  following  indentation.  The  Ni  tip 
occupies  the  topmost  eight  layers.  The  configurations  (a-f) 
correspond  to  the  stages  marked  0,  Q,  S,  U,  W,  and  X  in  Fig.  3b. 

Note  the  crystalline  structure  of  the  neck.  Successive  elongations 
of  the  neck,  upon  increased  separation  between  the  tip-holder 
assembly  and  the  substrate,  occur  via  structural  transformation 
resulting  in  successive  addition  of  layers  in  the  neck  accompanied 
by  narrowing  (i.e.,  reduction  in  cross-sectional  area  of  the  neck). 
Distance  in  units  of  X  and  Z,  with  X  =  1  and  Z  =  1  corresponding  to 
61.2  8. 
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Fig.  7.  von  Mises'  shear  stress  (vJT,)  corresponding  to  the  configuration 

marked  T  in  Fig.  3b  (i.e.,  just  before  the  structural  transformation 
resulting  in  the  configuration  (d)  in  Fig.  6).  The  proximal  inter¬ 
facial  layers  of  Ni  and  Au  are  marked  by  arrows.  The  maximum 
contours  (2.9GPa,  marked  a)  occur  on  the  periphery  of  the  neck  (X,Z) 

=  (±  0. 1,0.3).  The  increment  between  contours  is  0.2GPa.  The 
contours  marked  h,  i,  J  and  k  correspond  to  1.1,  0.9,-  0.7  and  0.5GPa, 
respectively.  Distance  along  X  and  Z  in  units  of  X  =  1  and  Z  =  1 
corresponding  to  61.2  8. 
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